{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# _*Pricing Bull Spreads*_ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Introduction\n",
    "<br>\n",
    "Suppose a <a href=\"http://www.theoptionsguide.com/bull-call-spread.aspx\">bull spread</a> with strike prices $K_1 < K_2$ and an underlying asset whose spot price at maturity $S_T$ follows a given random distribution.\n",
    "The corresponding payoff function is defined as:\n",
    "\n",
    "\n",
    "$$\\min\\{\\max\\{S_T - K_1, 0\\}, K_2 - K_1\\}$$\n",
    "\n",
    "\n",
    "\n",
    "In the following, a quantum algorithm based on amplitude estimation is used to estimate the expected payoff, i.e., the fair price before discounting, for the option:\n",
    "\n",
    "\n",
    "$$\\mathbb{E}\\left[ \\min\\{\\max\\{S_T - K_1, 0\\}, K_2 - K_1\\} \\right]$$\n",
    "\n",
    "\n",
    "as well as the corresponding $\\Delta$, i.e., the derivative of the option price with respect to the spot price, defined as:\n",
    "\n",
    "\n",
    "$$\n",
    "\\Delta = \\mathbb{P}\\left[K_1 \\leq S \\leq K_2\\right]\n",
    "$$\n",
    "\n",
    "\n",
    "The approximation of the objective function and a general introduction to option pricing and risk analysis on quantum computers are given in the following papers:\n",
    "\n",
    "- <a href=\"https://arxiv.org/abs/1806.06893\">Quantum Risk Analysis. Woerner, Egger. 2018.</a>\n",
    "- <a href=\"https://arxiv.org/abs/1905.02666\">Option Pricing using Quantum Computers. Stamatopoulos et al. 2019.</a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import numpy as np\n",
    "\n",
    "from qiskit import Aer\n",
    "from qiskit.aqua.algorithms import IterativeAmplitudeEstimation\n",
    "from qiskit.circuit.library import LogNormalDistribution, LinearAmplitudeFunction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Uncertainty Model\n",
    "\n",
    "We construct a circuit factory to load a log-normal random distribution into a quantum state.\n",
    "The distribution is truncated to a given interval $[\\text{low}, \\text{high}]$ and discretized using $2^n$ grid points, where $n$ denotes the number of qubits used.\n",
    "The unitary operator corresponding to the circuit factory implements the following: \n",
    "\n",
    "$$\\big|0\\rangle_{n} \\mapsto \\big|\\psi\\rangle_{n} = \\sum_{i=0}^{2^n-1} \\sqrt{p_i}\\big|i\\rangle_{n},$$\n",
    "\n",
    "where $p_i$ denote the probabilities corresponding to the truncated and discretized distribution and where $i$ is mapped to the right interval using the affine map:\n",
    "\n",
    "$$ \\{0, \\ldots, 2^n-1\\} \\ni i \\mapsto \\frac{\\text{high} - \\text{low}}{2^n - 1} * i + \\text{low} \\in [\\text{low}, \\text{high}].$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# number of qubits to represent the uncertainty\n",
    "num_uncertainty_qubits = 3\n",
    "\n",
    "# parameters for considered random distribution\n",
    "S = 2.0 # initial spot price\n",
    "vol = 0.4 # volatility of 40%\n",
    "r = 0.05 # annual interest rate of 4%\n",
    "T = 40 / 365 # 40 days to maturity\n",
    "\n",
    "# resulting parameters for log-normal distribution\n",
    "mu = ((r - 0.5 * vol**2) * T + np.log(S))\n",
    "sigma = vol * np.sqrt(T)\n",
    "mean = np.exp(mu + sigma**2/2)\n",
    "variance = (np.exp(sigma**2) - 1) * np.exp(2*mu + sigma**2)\n",
    "stddev = np.sqrt(variance)\n",
    "\n",
    "# lowest and highest value considered for the spot price; in between, an equidistant discretization is considered.\n",
    "low  = np.maximum(0, mean - 3*stddev)\n",
    "high = mean + 3*stddev\n",
    "\n",
    "# construct circuit factory for uncertainty model\n",
    "uncertainty_model = LogNormalDistribution(num_uncertainty_qubits, mu=mu, sigma=sigma**2, bounds=(low, high))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZwAAAEyCAYAAADOV2anAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3debgcVZnH8e9PEAgEMLKERSSACiJRMEGMoCSArDpBEKKiYxCNjAKOIgwiQgAdAWVxYBQjakRGoyKiIIhhSRDZ40JYomwhBhQEAjGEQELe+eOcC51K973dt/tWNcnv8zz93FtVp6re7tu33646myICMzOzgfaKqgMwM7OVgxOOmZmVwgnHzMxK4YRjZmalcMIxM7NSOOGYmVkpKk84kraVdI2khZIekXSKpFX62OdNkn6Tyz8naY6kCyRtXCg3WVLUeWwzsM/KzMyKVq3y5JKGAFcDdwNjga2AM0mJ8IRedl0XeBC4EHgE2AI4CRghaceIWFJTdhZwaGH/2c3Et/7668ewYcOaKTognnnmGdZaa63Kzt9It8YFjq0/ujUucGz9UXVcM2bMeDwiNqi7MSIqewBfAOYB69SsOxZYWLuuyWO9GwjgrTXrJgO39ze+ESNGRJWuu+66Ss/fSLfGFeHY+qNb44pwbP1RdVy9feZWfUttH+CqiJhfs24KMAjYtcVjPZF/rtaJwMzMrLOqTjjbkG55vSgi5pCucPqsZ5H0CkmrSdoaOA24Dbi1UGxbSfNzXc8NklpNZGZm1gGKCsdSk7QYOCYizimsnwtcGBHH97H/b4C98uIMYN+IeKxm+2eA50l1RBsARwMjgF0iopiYevaZAEwAGDp06IgpU6b056l1xIIFCxg8eHBl52+kW+MCx9Yf3RoXOLb+qDquMWPGzIiIkXU3NrrXVsYDWAz8Z531c4H/bmL/1wM7AR8mXSnNANbopfyapMYGlzYTn+tw6uvWuCIcW390a1wRjq0/qo6LLq7DmUdqcVY0JG/rVUTcGxG3RMRFpCudHYAP9VJ+IXAF8Nb+hWtmZv1VdcKZRaGuRtJmpCuRWXX3aCAiHgKeBLbsq2h+mJlZiapOOFcCe0lau2bdOOBZYHorB8oNB9Yj3TJrVGYQsB/p1puZmZWo0o6fwPnAUcAlkk4nXZ1MBM6KmqbSku4DpkfEYXn568AS4BbgKeCNpP4795OaVSNpXeBy4CLgPmB94LPAJsBBJTw3MzOrUWnCiYh5knYHzgMuIyWPs0lJp9aqQO1wN7cDR5Jak60BzAF+Dnw1Ip7JZZ4D/kkasWBDYBFwE7BrRNw+EM/HzMwaq/oKh4i4G9itjzLDCstTyFcyveyzCDig3fjMAIYd9+u2j3H08CWMb/M4s0/br+04zKpSdR2OmZmtJJxwzMysFE44ZmZWCiccMzMrhROOmZmVwgnHzMxK4YRjZmalcMIxM7NSOOGYmVkpnHDMzKwUTjhmZlYKJxwzMyuFE46ZmZXCCcfMzErhhGNmZqVwwjEzs1I44ZiZWSmccMzMrBSVJxxJ20q6RtJCSY9IOkXSKn3s8yZJv8nln5M0R9IFkjauU3aspJmSFkm6W9K4gXs2ZmbWyKpVnlzSEOBq4G5gLLAVcCYpEZ7Qy67rAg8CFwKPAFsAJwEjJO0YEUvy8XcBfg58EzgK2Bf4saR5EfHbAXlSZmZWV6UJBzgcGAQcEBHzgamS1gEmSjojr1tORNwI3FizapqkucBvgTcDf8jrvwRcHxFH5eXrJL0JODGXNTOzklR9S20f4KpCYplCSkK7tnisJ/LP1QAkrQ6MAX5aKDcFGCVp3dbDNTOz/qo64WwDzKpdERFzgIV5W68kvULSapK2Bk4DbgNuzZu3Al5ZPD5wD+l5v6G90M3MrBWKiOpOLi0GjomIcwrr5wIXRsTxfez/G2CvvDgD2DciHsvbdgZuAHaIiD/V7PM64F5gr3r1OJImABMAhg4dOmLKlCn9fXptW7BgAYMHD67s/I10a1wwcLHNfPjpto8xdBA8+mx7xxi+aecvzFfGv2cndGtsVcc1ZsyYGRExst62qutw2nUk8Grg9aRGBldK2jkiFvX3gBExCZgEMHLkyBg9enQn4uyXadOmUeX5G+nWuGDgYht/3K/bPsbRw5dw5sz2/uVmHzK67TiKVsa/Zyd0a2zdGhdUn3DmkVqcFQ3J23oVEffmX2+R9DtSy7UPAd+r2b94/CE15zYzs5JUXYczi0JdjaTNgDVZvu6lVxHxEPAksGVedT+wuHj8vLwU+Gs/4jUzs36qOuFcCewlae2adeOAZ4HprRwoNxxYj3SVQ0Q8B1wHHFQoOg64KSLavylvZmZNq/qW2vmkDpmXSDqddHUyETirtqm0pPuA6RFxWF7+OrAEuAV4CngjcCzpqqa2lv9UUh+dc4BLSR0/9wX2HtinZWZmRZUmnIiYJ2l34DzgMlLyOJuUdGqtCtQOd3M7qcHABGANYA5pRIGvRsQzNce/QdL7gS8D/0Gu4/EoA7aiGNahxgztNoqYfdp+bcdhK76qr3CIiLuB3fooM6ywPIVlr2R62/dS0tWNmZlVqOo6HDMzW0k44ZiZWSmccMzMrBROOGZmVgonHDMzK4UTjpmZlcIJx8zMSuGEY2ZmpXDCMTOzUjjhmJlZKZxwzMysFE44ZmZWCiccMzMrRcujRUsaDrwN2Ig0NcCTpNkzb4wIT9tsZmZ1NZVwJG1Jmk/mEGAoaYrmp4DngFeRpoReKmk6cAHwk4hYOiARm5nZy1Kft9QkXQDcBWwPnALsAKwRERtExGsiYjCwIfBeYCZwBnCPpF0GLmwzM3u5aeYK51lgm4h4qFGBiHgcuBK4UtLngIOATTsTopmZrQj6vMKJiCN7SzZ1yi+NiJ9ExE+aKS9pW0nXSFoo6RFJp0hapY99dpT0fUn35f3+IukkSWsUyk2UFHUeezf7fMzMrDPammJa0nbAroCA6RExs8X9hwBXA3cDY4GtgDNJifCEXnYdl8ueDtwLvBk4Nf88sFD2aaCYYO5pJU4zM2tfvxOOpP8AvgJcA6wFfE3S0RHxzRYOczgwCDggIuYDUyWtA0yUdEZeV89p+TZej2mSFgHflrR54YpsSUTc3EJMZmY2AJppNLBmg03/BYyKiIMiYl/gCOCLLZ5/H+CqQmKZQkpCuzbaqZBsevwx/9ykxRjMzKwEzXT8/KukQ+qsF6l5dI/+NIPeBphVuyIi5gAL87ZWjMox3F9Y/ypJj0taLOmPkg7oR5xmZtYmRUTvBaR3AecAzwNHRcStef2nSc2kryH1w9kdODYizm365NJi4JiIOKewfi5wYUQc3+RxNgLuAK6IiPE16z9MarL9R2Bt4JPAvsCBEXFJg2NNACYADB06dMSUKVOafTodt2DBAgYPHlzZ+Rvp1rhg4GKb+fDTbR9j6CB49Nn2jjF803WXWe7WuDplZXyvtavquMaMGTMjIkbW29ZnwgGQJOAwUsX8VOC/IuLvkt7CS7e+ro+IP7USWCcSjqTVSA0PXgOM6G20g/w8bgQGRcT2fR175MiRcfvtt/dVbMBMmzaN0aNHV3b+Rro1Lhi42IYd9+u2j3H08CWcObOtdjrMPm2/ZZa7Na5OWRnfa+2qOi5JDRNOU2OpRXIBsDXwKHCnpC8CsyLif/KjpWSTzQPqfTUakrf1KieQC4E3Afv2NbROpOx6CfDmvppem5lZZ7U0eGdEzI+IY4CdSOOpzZL0/jbOP4tCXY2kzUi36GbV3WNZ55CaU4+NiGbKA0R+mJlZiZpqpSbpy5JuyZXuk4BFETGWVNdxkqTp+fZaq64E9pK0ds26caTRDab3EdcXSC3jPhwRNzRzsnxFdCDw54h4oR/xmplZPzVz4/a7wLakPjcLSUlmqqRtI2JqTjSfyusujYgJLZz/fOAo4BJJpwNbAhOBs2qbSku6j9Sx9LC8/CHgv4HJwMOS3l5zzPsj4p+53HTg56SrpbWAT5CuzvZvIUYzM+uAZhLOPsBBETEVQNLvgSdIPf3vy6NCnyfpR8BJrZw8IuZJ2h04D7iMNAL12aSkU4yzts5lz/xzfH7UOpSUiADuA/4T2JjUZPoPwH4RcWUrcZqZWfuaSTizgI9ImgEsIjUtfgaYW1soIp4EPtNqABFxN7BbH2WGFZbHs3yiqbffYa3GY2ZmA6OZhPNR0hXD46TK9tmkK55FAxeWmZmtaPpMOBHxF2CUpLWA1Tyrp5mZ9UfTvb0i4hnSrTQzM7OWNdMs+iOtdpKU9DpJ7+x/WGZmtqJppuPn54D7JZ3aW18bSetJOkTSZcCfSC3DzMzMgObqcHaQNA44EviipAWkCcweB54DXgVsAbyWNBzNRcDhEfHwgEVtZmYvO03V4eTpon8iaStgD+CtwEakzpSPAtcDvwemRcTiAYrVzMxexloaIjYi7mf5+WbMzMz61NLgnWZmZv3lhGNmZqVwwjEzs1I44ZiZWSlaSjiS3ivJScrMzFrWavK4FJgr6XRJbxyIgMzMbMXUasLZCvgOcDBwp6SbJH1C0jqdD83MzFYkLSWciJgdESdFxBbAu0kTnJ0N/F3SDyWNGYggzczs5a/f9TERcW1EfAR4AzADOAS4WtIDkj4rqaVOpWZmtmLrd8KRtKukycBfgO2A/yVN/XwxcDJwYScCNDOzFUOrrdQ2l3SipPuBa4HNgAnAxhFxZERcExHHkmYJHdvkMbeVdI2khZIekXRKX9MhSNpR0vcl3Zf3+4ukkyStUafszpJukbRI0oOSjmrlOZuZWWe0etvrAeAR0pTT34uIBxuUuwu4ta+DSRoCXA3cTUpQWwFnkhLhCb3sOi6XPR24F3gzcGr+eWDN8V8HXAVcDnwBeBtwlqSFEXFBX/GZmVnntJpw3gNcFRFLeysUEX8FmmlAcDgwCDggIuYDU3OLt4mSzsjr6jktIh6vWZ4maRHwbUmbR8RDef0xpAT54YhYAlwr6bXASZK+GxHRRIxmZtYBrdbh7EialmA5kjaWdGKLx9uHlMBqE8sUUhLatdFOhWTT44/55yaF41+Sk03t8V9DqncyM7OStJpwTiJ9WNezSd7eim2AWbUrImIOsDBva8UoYCl5+gRJa5HqmGYVyt1Tc24zMyuJWrmrJGkpsFNE3FZn21jguxGxfgvHWwwcExHnFNbPBS6MiOObPM5GwB3AFRExPq/bFJgLvC8iLq0puyqwGPhkREyqc6wJpIYQDB06dMSUKVOafTodt2DBAgYPHlzZ+Rvp1rhg4GKb+fDTbR9j6CB49Nn2jjF803WXWe7WuDplZXyvtavquMaMGTMjIkbW29ZnHY6kj5JanQEE8C1JxbqVNYDhwG/bCbQ/JK0G/BRYAHy23ePlJDQJYOTIkTF69Oh2D9lv06ZNo8rzN9KtccHAxTb+uF+3fYyjhy/hzJntdU+bfcjoZZa7Na5OWRnfa+3q1riguUYDC4En8u8CngaeLJR5HrgS+GaL558H1PtqNCRv65Ukkfr7vAnYOSJq93kq/ywef0jNuc3MrCR9JpyI+BnwMwBJ3wdO6aU5dKtmUahLkbQZsCbL173Ucw6pOfW7I6JYF/SMpL8Vj1+z3MzxzcysQ1odS+3QDiYbSFdFe0lau2bdOOBZYHpvO0r6AnAEqcnzDb0c/32FjqTjgL8Bd/Y7ajMza1nVc9ucDzwHXCJpj1xhPxE4q7apdB5R4Ls1yx8C/pt0O+1hSW+veWxQc/yvkVrV/VDSGEnHAp8kXaW5D46ZWYmaaTRwKzA+Iu6WdBup4UBDEfG2Zk8eEfMk7Q6cB1xGqnc5m5R0inHWXqXsmX+Oz49ah5JGQiAi7pO0N3AW6WrnH8DRHmXAzKx8zTQauIt0i6vn945eGUTE3cBufZQZVlgez/KJptG+N5CGtDEzswo102jg0Jrfxw9oNGZmtsKqug7HzMxWEs3U4fRZb1OrlTocMzNbeTRbh+MWXWZm1pZm6nDGlxCHmZmt4FyHY2Zmpai0H46Zma08Ku+HY2ZmKwf3wzEzs1K0PAlGnn9mPKn3/sbA34FbgB9ExPMdjc7MzFYYLTUakPRG4F7gf4HtgBfyz/8F7pO0bccjNDOzFUKrVziTSBOwvTMi5vSslPRa4HLS6M/v6lx4Zma2omg14YwEPlibbAAiYo6kk4AfdSwyW+kM69B0ye1Ouzz7tP3ajsPMltdqP5zZwBoNtq0BzGmwzczMVnKtJpzjgC9L2ql2paS3A6cC/9WpwMzMbMXSn8E71wFulPQY8BiwYX48ARwPXDoAcZqZ2ctcfwbvvGuAYjEzsxVY5YN35qbU5wKjSFNMXwCcHBEv9LLPasBXgLeTGjKsERGqU24y8NE6h3hjRMxqP3ozM2tWyx0/O0nSEOBq4G5gLLAVcCapbumEXnZdE/g4cCtwI71PUT0LOLSwbnb/IjYzs/6qNOEAhwODgAMiYj4wVdI6wERJZ+R1y4mIpyS9OiJC0hH0nnCeiYibOx+6mZm1ouXpCSSNk3S1pDmSHis+WjzcPsBVhcQyhZSEdu1tx4jwIKJmZi8jrQ5t8yHgB8B9wGuAX5FGGHgFMB84r8Xzb0O65fWi3Kl0Yd7WCdtKmi/pOUk3SOo1kZmZ2cBQKxcKkv4IXAycBiwGRkbEHyStDUwFLo6Ir7dwvMXAMRFxTmH9XODCiDi+iWMcAZzboNHAZ4DnSXVEGwBHAyOAXSLi1gbHmwBMABg6dOiIKVOmNPt0Om7BggUMHjy4svM3MlBxzXz46baPMXQQPPps3+V6M3zTdZdb162xdWtcndKt/wPQvbFVHdeYMWNmRMTIettarcN5PfD7iHhB0gukPjlExL8knQ6cDTSdcAZaRHyjdlnSFaRm3ccD+zfYZxJpzDhGjhwZo0ePHuAoG5s2bRpVnr+RgYqr3SFpIA1tc+bM9qomZx8yerl13Rpbt8bVKd36PwDdG1u3xgWt1+HMB1bPvz8MvLFmm4D1WjzePKDeV6MheVtHRcRC4ArgrZ0+tpmZ9a7VrzW3AW8GriLV35woaQnpttWJQKutwWZRqKuRtBmp2fNA9ZMJPGupmVnpWk04XwU2z7+fmH//FulK6Tbgky0e70rgGElrR8S/8rpxpCmtp7d4rD5JGgTsB8zo9LHNzKx3LSWc3J/l5vz7U8BYSasDqzfqM9OH84GjgEtyHdCWwETgrNrjSboPmB4Rh9Ws2wdYC9g+L78/b7otIh6StC6pBd1FpFZ16wOfBTYBDupHrGZm1oaOTTEtqeUppiNinqTdSc2pLyMNbXM2KekU41ylsO5bvHS1BfCz/PNQYDLwHPBP0ogFGwKLgJuAXSPi9lbiNDOz9rWUcPIU078hXSXMII0WvR3w78CXJO0dEXe3csxcvreRAoiIYc2sK2xfBBzQSixmZjZwPMW0mZmVotVm0SOBE+tNMQ2cBOzYqcDMzGzF4immzcysFK3eUjsOOFPSgxFxS8/KmimmP9/J4Mzs5WtYh0ZBaHc0hdmn7dd2HNYZnmLazMxK4SmmzcysFJVPMW1mZiuHfg0RK2kTYBTwatKttJsj4pFOBmZmZiuWVjt+rgKcC3yCZXv+vyBpEnBkRCztYHxmZraCaLVZ9MnAx0iNA4aRpoIelpc/xvJD0piZmQGt31L7d+CEwqyec4CvSQrSQJwndio4MzNbcbR6hbMhcEeDbXfk7WZmZstpNeH8FfhAg20fAP7SXjhmZraiavWW2peBKXmwzouBR0lXNQcBY2icjMzMbCXX6gRsP5X0FKnxwDeAVwKLSVMV7B0RUzsfopmZrQiaTjiSXkmadO3OiBgl6RWkWTQfd1NoMzPrSyt1OC8A1wLbAETE0oh4zMnGzMya0XTCyYnlXmCjgQvHzMxWVK22UvsicKKk4Z0KQNK2kq6RtFDSI5JOySMa9LbPapK+Jul3kp7NfYAalR0raaakRZLuljSuU7GbmVnzWm2ldgKwHvAnSQ+TWqkt82EfEW9r9mCShgBXA3cDY4GtgDNJifCEXnZdE/g4cCtwI7Bbg+PvAvwc+CapU+q+wI8lzYuI3zYbp5mZta/VhHMXcGcHz384aXicAyJiPjBV0jrAREln5HXLiYinJL06IkLSETRIOMCXgOsj4qi8fJ2kN5FGQ3DCMTMrUavNosd3+Pz7AFcVEssU4HRgV+CyXmJpeBsNQNLqpL5BRxU2TQG+L2ndiHi6X1GbmVnLmqrDkTRI0oGSjpb0IUlDO3T+bYBZtSsiYg6wMG9rx1akfkKzCuvvIT3vN7R5fDMza4H6uFBA0pakepZhNavnAwe3Ww8iaTFwTEScU1g/F7gwIo5v4hhHAOdGhArrdwZuAHaIiD/VrH8dqbXdXvXilzQBmAAwdOjQEVOmTGn9iXXIggULGDx4cGXnb2Sg4pr5cPsXnEMHwaPPtneM4Zuuu9y6bo2tW+OC7o6tE1a2/89mjRkzZkZEjKy3rZlbamcAS4F3kkYU2IJUCf/t/PsKJSImAZMARo4cGaNHj64slmnTplHl+RsZqLjGH/frto9x9PAlnDmzX/MKvmj2IaOXW9etsXVrXNDdsXXCyvb/2QnN3FIbRZqS4PcRsSgi7gE+CbxW0sZtnn8eUO/rx5C8rd1jU+f4QwrbzcysBM0knI2BBwrr7gdE+51AZ1Goq5G0GanZc7HupVX3k8Z5K9YFbUO6Yvtrm8c3M7MWNNvxs/eKnv67EthL0to168YBzwLT2zlwRDwHXEcaybrWOOAmt1AzMytXszdHr5K0pM76a4rrI6KVSdjOJzVbvkTS6cCWpGmqz6ptKi3pPmB6RBxWs24fYC1g+7z8/rzptoh4KP9+KjBN0jnApaSOn/sCe7cQo5mZdUAzCefkgTp5RMyTtDtwHqnPzVPA2aSkU2tVoDjczbeAzWuWf5Z/HgpMzse/ISeiLwP/ATwIfMijDJiZla/PhBMRA5Zw8vHvpvFIAT1lhjWzrsG+l5KubszMrEKtDt5pZmbWL044ZmZWCiccMzMrhROOmZmVwgnHzMxK4YRjZmalcMIxM7NSOOGYmVkpnHDMzKwUTjhmZlYKJxwzMyuFE46ZmZXCCcfMzErhhGNmZqVwwjEzs1I44ZiZWSmccMzMrBSVJxxJ20q6RtJCSY9IOkVScTrpevutK+n7kuZJelrS/0lar1BmsqSo89hm4J6RmZnV0+cU0wNJ0hDgauBuYCywFXAmKRGe0MfuPwXeAHwcWAqcTppK+p2FcrOAQwvrZrcTt5mZta7ShAMcDgwCDoiI+cBUSesAEyWdkdctR9IoYE9g14i4Pq97GLhF0h4RcXVN8Wci4uaBfRpmZtaXqm+p7QNcVUgsU0hJaNc+9nu0J9kARMStwIN5m5mZdZmqE842pFteL4qIOcDCvK3p/bJ76uy3raT5kp6TdIOk3hKZmZkNEEVEdSeXFgPHRMQ5hfVzgQsj4vgG+00l3Srbv7D+ImDLiHhHXv4M8DypjmgD4GhgBLBLviKqd+wJwASAoUOHjpgyZUobz7A9CxYsYPDgwZWdv5GBimvmw0+3fYyhg+DRZ9s7xvBN111uXbfG1q1xQXfH1gkr2/9ns8aMGTMjIkbW21Z1Hc6Aiohv1C5LugK4Czge2L/BPpOASQAjR46M0aNHD3CUjU2bNo0qz9/IQMU1/rhft32Mo4cv4cyZ7b2tZx8yerl13Rpbt8YF3R1bJ6xs/5+dUPUttXlAva8fQ/K2ju4XEQuBK4C3thCjmZl1QNUJZxaFOhdJmwFrUr+OpuF+WaO6nVqRH2ZmVqKqE86VwF6S1q5ZNw54Fpjex34bSdqlZ4WkkcCWeVtdkgYB+wEz2gnazMxaV3XCOR94DrhE0h65wn4icFZtU2lJ90n6bs9yRNwE/Ba4UNIBkvYH/g+4oacPTh6J4HeSPilpd0njgOuATYD/LusJmplZUmmjgYiYJ2l34DzgMuAp4GxS0qm1KlAc7mZcLvs9UuK8HDiqZvtzwD9JIxZsCCwCbiJ1Fr29o0/EzMz6VHkrtYi4G9itjzLD6qx7ijRkTXHYmp7ti4ADOhCima1ghnWoBV07LfFmn7Zf2zG83FR9S83MzFYSTjhmZlYKJxwzMyuFE46ZmZXCCcfMzErhhGNmZqVwwjEzs1I44ZiZWSmccMzMrBSVjzRg5eqGHtawcvayNlvZ+QrHzMxK4YRjZmalcMIxM7NSOOGYmVkpnHDMzKwUTjhmZlYKJxwzMyuFE46ZmZWi8o6fkrYFzgVGAU8BFwAnR8QLfey3LnAOsD8pcV4OHBURTxTKjQW+DLweeCAf+yedfh5mZu1a0TtmV3qFI2kIcDUQwFjgFOBo4OQmdv8pMBr4ODAe2BG4tHD8XYCfA9cB+wC/Bn4sac+OPAEzM2ta1Vc4hwODgAMiYj4wVdI6wERJZ+R1y5E0CtgT2DUirs/rHgZukbRHRFydi34JuD4ijsrL10l6E3Ai8NuBe1pmZlZUdR3OPsBVhcQyhZSEdu1jv0d7kg1ARNwKPJi3IWl1YAzpSqjWFGBUviVnZmYlqTrhbAPMql0REXOAhXlb0/tl99TstxXwyjrl7iE97zf0I14zM+snRUR1J5cWA8dExDmF9XOBCyPi+Ab7TQWeiYj9C+svAraMiHdI2hm4AdghIv5UU+Z1wL3AXhGx3G01SROACXlxa+Av/X6C7VsfeLzC8zfSrXGBY+uPbo0LHFt/VB3X5hGxQb0NVdfhdJ2ImARMqjoOAEm3R8TIquMo6ta4wLH1R7fGBY6tP7o1Lqj+lto8oF5dypC8rZ39en4Wyw0pbDczsxJUnXBmUairkbQZsCb162ga7pfV1u3cDyyuU24bYCnw137Ea2Zm/VR1wrkS2EvS2jXrxgHPAtP72G+j3M8GAEkjgS3zNiLiOVL/m4MK+44DboqIp9sPf8B1xa29Oro1LnBs/dGtcYFj649ujavyRgNDgLuBO4HTSQnjLOCciDihptx9wPSIOKxm3VWk0QM+T7piOR14LCLeWVNmF2AacB6pU+i+ufze9RoMmJnZwKn0Cici5gG7A6sAl5FGGDgbOKlQdNVcptY40lXQ94ALgRnA+wrHvwF4P7AHcBXwb8CHnGzMzMpX6RWOmZmtPKquwzEzs5WEE1I1XdsAABkPSURBVI6ZmZXCCcfMzErhkQa6gCSRGjzsB7wReDWp5d0/gJuByRFRSb+h3C9qX0DAzyLiCUmvIbX22wqYDUyKiJklxvRfwBVlnrNZkgYBq0bEv2rWbQAcAWxL+rv+Cfjmy6RpfiXy/8R7gbeSpi+5nfQ374pK5zyq/ePAbrlxUlUx7AasBvw6Ip7J77VPk1r8PkD633ykivjqcaOBiuU3yBXACFKCeR7YlPRPdiXpjbM1cGpEnFpybG8DpgJrAUuAJ4G9crwvAHcB2wEbAXtExO9Kimsp6fWZBfwImBIR95dx7r5IugK4NyI+k5dHkf6OS0ktKUX6Wz9P+rC6q6S4dgAGRcSNNev2Br7AS4nwz8DE2jIlxXYjcFhE3JOXh5CmDxkBLMjFBpO+fO1Vm8wHOK5P9bJ5EPA14BuksRmJiG+WERe8OCbkNcBmedWDpClbpgKvInV835rUp3FERMwtK7ZeRYQfFT6AH5PesMNr1m0C/Ab4eV7elfSP97GSY5tK6jz7KtLI2+cBc4FfAq/MZVYnfaBeV2JcS4GvkmZ5fY6U/G4DPgtsWvHf83FgbM3yzaQPhrVr1q1LatJ/VYlx3Qx8sWb5Y/l1vAb4InBC/lsvqY2/xL/n22qWv0v6crN3zbq9ScNRnV1yXC/kn/UetdteKPk1+ynpC8LrSHdEfpg/R27sea+RBvH8M/DtMmPrNe6qA1jZH6RptQ+ss35YfkNvnJePB/5ccmxPAPvULG+Y/7n2LJTbD3i8xLhe/IDK/2wT8gfnkvyYlte9uoK/50LgXTXLzxdfr5rX7JkS45pfGwdwH3BunXLnV/A+KyacfwL/Wafc54GHSozrF8DfgUPJd4Nqtr0qx/2usuIpnP8R4OCa5c1zPAcUyh0K/LWKGOs93Gigeq8gJZaiF0i3X3oGH72FaubwiTq/F+/DVnZfNiKejIhJEbE78BrSFOWrkT44/y6p/UniW3MnaeK/Ho+SkmLReqTkVJalheXNgYvrlLuYdCumSq8i1dkUzSDdvi1FRLwP+ChwDHBbnvLkxc1lxdHAENIt+B4P558PFco9QPq/6ApOONWbCnxZ0pY9K/I97P8hvaF6GgsMBsquZL4d+LykwZJeQbrKehj4D0mr5FhXBT5F+qCtVET8IyK+ERHvALYgjVixSclhnAYcJ+lj+bX5CvA1Se+WtJqk1XPdyVdZfjbagfQ74JCa5buAekPY78hLH15lOlDSp3K9yTyg3nwq65Ou1EoTaVSSN5NmCv61pCm50UzVHiN9aejxAvBt0hecWhsCpdR5NaXqS6yV/UH69nEXaWTr+0hjyz1LutVWezvrDOAnJcc2kvTPvzjH9ATwFlISfIA0HNGDpHqUMSXGtcwtmG57AB8nfTA+Ddyaf3+BdLvvhfz4BbBmiTENz3H8EHgbaSr2x0gJ8d2kCufTgEXUuZ1Vwt+z+PhenXLfBn5X4d91I9IwWguAM/PfsapbapfWe43qlDsXmFrVa1Z8uJVaF8hXCweTPszXICWeH0XEk5UGBuRvc+8hNaH/eUT8XdJGwLGkWy8PARdExB9KjOkk4DvRRc09iyStRxrv722kDyqRkvc9wOURMaOCmLYHvgXsRLolpLyp5/d5wCkR8Y2yY2uGpE8A90fEtRXHMYo05uPWwH5Rcqu+HMNQ0heWB/so9zlSndw15UTWOyccs5WMpDeSkk4xEd4YEYurjM1WbE44XUTSm0gTxNXOSjorSuqr0SpJr4iIYmV0ZSStQeqMuhS4r+oPz1yHsyU1HXkjYk6VMb3c5A6gRIUfVLkzryJiYc267ckdn6u4Wn25cqOBLpArmB8C7gB+RppAaVL+/Q5JsyUdWlFsB0i6VNIVkt6b142TNBtYLOmhfKujzJg+LOljNcurSjqN1HfjDlIDhiclHVdmXDXxjJD0K1Jl7T3A74GbgAclPSzpFElrVhFbN5K0Z2ESRiTtL+kPpPrD5yXdLmm/kuNaV9IvSHVf8yV9R9Iqkn4A/IH0/3mrpBskrV9mbM2SdKCkeq1gK+GEUzFJR5IqQy8HRpNalbwyPzYkdfq8HDhf0qdLju1gUjPZ9Un/+D/JyeWHpH4vR5E6mp0vaa8SQzue1OG0x+k5lq8C7yK9ZmcCJ0k6vsS4kLQn6TXZhHSf/1RSr/kXgImkCQYPBG7MrRHLjO09kq6RdI+kX0p6V50yO1XwAXUlaUinnhjeB1xCasBwXH4sBn6ZX9+ynAq8E/gcqaPsO0gtC3cjdUQdSqrf3CKXtT74llrFJD0AnB8RZ/RR7ljg8IjYsrdynSTpNmBGRByelw8hTXh3XkQcXVPu+8BmEbFHSXEtJLXgm56XHwO+UqzslvR54MiI2LzOYQYqthnAnRHx0cL6I0l9hLYk9RO6Ebg5InobPqWTcb2bNHrFzcAfgVHA9sA5wOd7bllJ2olUl1Oc8HAgY1sKvD0ibs3LfwAejoj3FspdAawVEbuWFNds0vvqO3l5B1JfoEMj4gc15T4BHB8RW5QRVz7n95osujkwusy/Z298hVO9jUlNZ/tyKyV2esu2ZtnOgZeTrryKnSkvIY3HVZanSVddPdYlDeFR9GfSVWKZtgUuqrP+IuC1wNYRsYj0Qf++OuUGyknAhRGxc0QcEREjgE8AnwQuyfVf3WI70lV/0STSYJ5l2YCX+sFBHjONNE5Zrfuo329oIH2U1JR9eB+P0r5sNcMJp3p/Bj6RO1bWlStOP0GqnyhTsOzU3j0DKT5VKLeA1Du8LL8idUhdLS9fDXywTrkPkvo4lekxUvP2oreQXs+ezrsP8dIoEmXYjkIijIjvkW4/vh24VlK9ERHKUnur5Wnqd1Z8hnI/sx4kvT493klq/PGOQrmdgbIbg9wLXBsRO/b2IN2O7BqenqB6R5Nuddwt6RLSCMg9H+jrklqtvY/UQXTvkmN7iPSN/SqAiHgh90G4p1BuK9KYU2X5Aqnn/J2SLiB1QD1d0nakcdREGl5mB9IQ92WaBJwqaTApET5P6r3/RdIApz19h7ak3A+pRaRRv5cRETPykC1XkW7zTSwxplpXSVqSf1+X9LebXiizDeW+z84HviFpOCkJHkx6730p/33/TLri+izl1+HczPKJr57a/laVc8KpWET8PjexPJY09MhmhSJ/I1Wqfi3KH4L/EgrjMEXELXXKfQAobU6QiHhS0ttJH+Kf46XbZqPy43nSkEHvjIjbyoorx/aVXCdxHHBiz2rSqOD/WVN0MfDfJYZ2B2l0gV8VN0TEAznpXAFMLjGmHifXWfdYnXUHkka0LkVEnJfvPHyQ1DDg2Ig4X9Jc0tBTPePhnQ98vay4snNJLeX6Mp1lx/arlBsNdJncXLbn9tRTtW3/u5Wk15JiLXWcq5rzD2PZToz3d0EfnFeSrvzWAB6o6rWpieeTpNZ9OzQawULSWqQhd/aICN9u70W+zb1+RPyz6lheTpxwzMysFL6l1iWUpnLehPTt/PE629cH9o2IC0sPro58D/sPwCFl37ZSl0/jrC6clrvb6WUyXXK+sqmd+noGKd7Sv7lLGkm6zSjSNPSzJL2FdIuy5312XkRcVXZsjfgKp2KSVie1Hjogr1pKGpH2c7UflhX1j9i3l81rAT8h1VXcCRARV5QUV1dO45xj6cppuZulNM7aQRFxSonn7MrpktWlU1/nWPYiNZZ5ktR6bwNgLKne9R5SX6sRpAYrB0bEpWXF1qsyhqT2o9fhw08ktUr7BGk6gKNIc1rcC7y+ptxOlD+NbVdOsUuXTuOcz9uV03K3EP+BFbzPunK6ZLp06ut83t+ThtZZJS8fn+P4bqHcD0kdjCt/b0V4iunKH6Rm0EcU1m0EXE+aandUXldFwrmdl6bY3bzweHP+hzy4Z12JcXXlNM75nN06Lfdrm3wcXsH7rCunS66TcLpi6ut8zqdJV8g9y0NyvLsVyu1JatBTWmy9PVyHU73NKHTojIh/SNqd9O3k6jykTJn9D3rsSLryOp3U7+XzkeffkNTTafEfEVGc1nag9UzjfH1e7pZpnHt047Tcs5s8p5os10kvl+mSu2Lq6+xZlu1X1fP7oEK5NUl9sLqCE071HgFez0sfngBEatb7AUnfIF06l95YINJXpEmSfgp8GZgp6dz8e5VOA/5P0t9Ir0vPNM5PkG6j9XT8LHsaZ3hpWu4bSMmudlruayN1nq1iWu5/AdcCF/RRbhdSn7AydfN0yQfmynnooqmvSbfUTpR0bz7310mzBf+XpOsj4l/5S+GxpITYHaq+xFrZH6TBMKf1UeYLlFxP0iCON5N68j8CfIZqp9jtummcc1zdOi33VNJQKH2Vq6IOpyunS6aLp74m1XfNrnmv30+6JdrzvzCTlJznAduXGVtvD1/hVO+bwDhJr44GHfIi4qtK8+W8u9zQlovjDmC0pA8AZ1DhkBkRcUGeq6RnGucn6YJpnCPi9jwUSnFa7nfx0rTcV1LytNykK+gJTZT7J4Wr7RJ8knTrpy8PkpJTKaL5zq+3k1psliYi7stDOe1MapxyTUQ8K2k06cvY1qRb8j+Kklr1NcPNoq1f8m2htYAFEdE1EzyZWfdywjEzs1J4vKSXiTy97XerjqOebo2tW+OC7o6tW0m6WtI1VcdR1K1xQffF5jqcl48xdO8XhG6NrVvjgi6NTdLVpDsfu1cdSx2iC18zujcu6LLYfEvNzF6Uvw2/IiK6Zkh7W3F0Teaz3klaI08D0HW6NbZujQu6N7aI2L1bk42kV3bja9atcUH3xeaE8/KxH6lZaDfq1ti6NS7o0tiq+oCS9GlJ90t6VtKfJX2kTrG3UvJr1q1xdXtsjTjhmK0kuvUDKvfrOpc0COuXSJ0YJ0u6WNIaZcbycoir22PrjRsNVEzStU0WrTekxoDq1ti6NS7o3thqPqB+TBq6/h2kD6ixwIcjosrxtj4PfD0iXhxSJ48l+H/AdZLeExFPOK6XTWwNudFAxSQtAf5CGgepN5sCO0W58+F0ZWzdGhd0b2ySbicNbVPvA+pB4D2RJoqrYt6lfwHvjYhphfXDSKMyrEKaBmCDMmPr1ri6Pbbe+AqnencBsyJiXG+FJL2fkofPoHtj69a4oHtj25r0rfhFEXGNpLeTPqBukrR3ifHUepo0AOYyImK2pHcAvwZuAk51XC/q5tgach1O9W4G3t5EuaD8scu6NbZujQu6N7aGH1Ck22uPkz6gdiwxph4zgP3rbYiIecDupPHK/qfMoOjeuKC7Y2vICad6ZwBHNlHuCmCLAY6lqFtj69a4oHtj6+YPqIuALSXVm9OIiHgW+DfS1ApzHBfQ3bE15Docs5WApIOAz5LqauqOSi5pFeBbwLsjouxEbSsBJxwzMyuFb6mZmVkpnHDMzKwUTjhmZlYKJxwzMyuFE471StJ4STMk/UvSPEl/lHTWAJ3rYEnjmyg3UVLUPB6R9HNJWzV5nsm5533lmn3OuWzP8763wfZ78/aJAxVDi8dd5nXu9HkkvULSEfk9+ayk+ZLukvQ/kvrVx0nJnyR9tMH2ybk3f71t58mT6vXKCccakvQFUjv+q4ADgH8Hfklq3z8QDgbGN1n2aWBUfnwe2B64RtJaTex7agvnGWitPGeARcAWkkbWrpS0IzAsbx/oGJpVfJ07fZ6fAF8GLiG9Jz9K6t/0juh/89uDgVcDP+rHvl8HDpH0un6ee4XnoW2sN0cA346I42vWXSbp5KoCqrEkIm7Ov98saQ7wO2Bf4GfFwrmPySoR8XxE3F9inJ32DPAH4AOkjpo9PgBcC4yoIqgeZb3OkvYB3g/sGxFX1mz6RX+vbrKjgB9GxOKac61KSp4fATYBPijpfuDkiHhxeKI8rMwNwH8AR7cRwwrLVzjWm1cB/yiurP322HPbRNL+kmZJWiTpBknbFvfLt1RmSnpO0t8kfSX/MyNpMnAgsGvNrbKJLcQ6I/8cVieuu0jf/Heq3VaI7V2SrpO0QNLTkqZJ2qFm+zslTZe0UNITkr4jae3eApI0StKvJP1d0jP5Vs0hta9dP5/zFODgng/W/PPgvL5jMeTX4OLC8UbnMtvVvpZ9vc6NziNpX0lLJW1ROM8Wef3YBq/BrvnncqNz9/fqJl+ZvAO4uLDpM8CxpFEYrgA+BnwPWK/OYX5OusrxZ2sdvsKx3vwBODJfPVzey3DnmwNnkebleBY4GbhK0ut7hr2XtCfpFsiFwDHAm0nfGtcDDs+/v5aU5D6Vjzu3hViH5Z//KKw7Azglr687z4uk0cBU4DrSbZlngJ1JIzr/UdLOwNXApaRv1esBpwFD8nIjmwO/B84nfRDvDHxf0tKI+DH9f86XkEYE2IV0VfdO0qjAlwBfKymGWsPo+3VudJ6/A4+QXveJNeXHA4+RBqGs55n882uSzoyIh1qMuZ7d83H/XFi/K2mk7TPyF6nf5zHo6rkRGAoMr3Mciwg//Kj7ICWFB0gDTS4ljYR8CrBOTZnJefs7atZtDiwBDq9ZdzNwXeH4xwIvAK/JyxcD05qIayJpsMlV8+MNpGQxH9i4ENf2dfafDNxes3wT6faUGpzvd3Vi3y0ff7smX0vlWL9N+vDqWd/Uc6593vn3XwL/m3//JnBp/v1xYGInYgCmARcX1o2ufd4tvs6NzvNlUpJSTZyzSfO9NHotNgLuyOcO4E7geGBwG+/3ScBtddZ/G/hbPudkYFgvx1g1v/c/0d84VuSHL/usoYi4A3gjqUL2m6QPgi8Bt0saXFP0sYi4sWa/h0i3uN4GL97XfyvL1638hHRbd1Q/wlsPWJwffwG2BMZFxN9ryjwcEX/q7SC5kcFOwA8if2IUtq+Z4/uppFV7HsAN+dwN60wkDVFqMfVQTawTSAmyXVOA90tanXSVtdzttBJi6NHn69yH75G+pIzOy2Py8vcb7RAR/wB2APYiXe29CvgKcKOk1QAkfTTfQvxTvo07K/8+Q9Ir6xx2I1LCLvoK6crnQdL/wufzVW+9uJYAT+VjWYETjvUqIp6LiMsi4oiI2Bb4OPB64LCaYo/V2fUxYOP8+/rAK4FHC2V6luuOeNuHp0lD6Y8EXkP61nlloUzxfPUMISXSv/eyfRVSwl1c83iO9Jw26+XYk4FxpNtce+Z4vwd0YgrgXwGDSR+GawGXVRBDj2Ze54Yi4gHS1dShedWhwK0RcVcf+70QEb+NiE+Rbtd9n3Qra1Te/oOI2J70ZWcJsHNEbB8RI6KmUUCNNUh/1+J55uTjvo90xb8LcIMadw94js6+visM1+FYSyLiu5LOALapWb1hnaIbkm7BQfrWuLhOuaH5Z93Ri/uwJCL66kvTTOXxPNLtwo0bbH8qH2ciqcK46JF6OynNK/8e4NMRcX7N+o58yYuIZyRdThoB+mcR8UyxTAdiWASsVlg3pF44TR6vNxcA31Fqin8ALbbyioilkn5LSlbFD/vXA/Oi7ymXn6TBlUlOUL9Rmqp7Immqh7MlnZMTUq1X0b/39ArPVzjWkKTlEomkDYB1WfZb7YZKswz2lHkt6VvlrZC+iZJusR1UONzBpA/7m/Ly85T8zTB/UN8C/HtPq686228Gto6I2+s86iYcYHXS/9eL35hzq7ZiH6Z2nvO3SFc25zfY3m4Mc1n2iwWkq6T+6u25XpK3TyHFXPcWIYCkoQ02/RuwkPT3rPUWmqvA/wt15iiq974Abss/X10ouwGwJvDXJs630vEVjvVmpqRfAr8l3SLbnNTJciHwg5pyjwMXSTqBl1qpPUa6ndPjJFLLte+TPkyGk1oufScielpFzQLGStqf9GH3SC8f6J10HKkV2pWSJpHu148iVXhfTmrccI2kpaSK73+RbuHsB3wxIpb7cImIpyXdBpwoaT4psR5HuhW4Tk3Rfj/nSPPZT+tle7sx/AI4TNLZpNZiY4B2pqFu+FwjYpGk/wM+Dfw4Ip7q5Tg/lfQv4KekxgUbAocAY0mV9cV930JqYNCX35Neqw0i4p81638k6Y/A9aTblyNIV5YPA/cUjjGSdMV3I7a8qlst+NG9D9I//29Jt40Wkf65fwRsU1NmMqmF1wGkb3XPkf5xl2u9RapLmEn6JjuXVP+was329Ukfck+Sb2M1iGsiubVWL7FPpqaFVF/bSE1frycl06dIrd62r9m+E/AbUku4Z4C7SU3B1+0lhtcB1+Tyc0iJa5nYm33OLTzvZVqptRsD8AVSC61/kWaZ/DeWb6XW1Ovc13MF9sjr9+jjOX4s/y3m5vfSk6SEOLpB+cuADzTxfl8NeAL4SGH9+/L5/kFK2vNJiX6HOsf4BoUWjX689PAEbNaW3KFvu4gY2VdZs97kusGDgS0jYmkHjzsH2Csiilcj9cp+A3hdROzXYPtkUqKcXWfbKsBDwHERcVFbQa+gfEvNzColaWtgW9KQMCd3ONkMIXWKbbZO5WvAXyW9IercKu3DQaRbyg3rn1Z2bjRgZlX7NulW7RWk4WM6JiLmRcSgSA1Xmik/l3TLrlGrxUtJt1zrEXBYpL44VodvqZmZWSl8hWNmZqVwwjEzs1I44ZiZWSmccMzMrBROOGZmVgonHDMzK4UTjpmZleL/AW46+1xJ9j1lAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot probability distribution\n",
    "x = uncertainty_model.values\n",
    "y = uncertainty_model.probabilities\n",
    "plt.bar(x, y, width=0.2)\n",
    "plt.xticks(x, size=15, rotation=90)\n",
    "plt.yticks(size=15)\n",
    "plt.grid()\n",
    "plt.xlabel('Spot Price at Maturity $S_T$ (\\$)', size=15)\n",
    "plt.ylabel('Probability ($\\%$)', size=15)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Payoff Function\n",
    "\n",
    "The payoff function equals zero as long as the spot price at maturity $S_T$ is less than the strike price $K_1$, then increases linearly, and is bounded by $K_2$.\n",
    "The implementation uses two comparators, that flip an ancilla qubit each from $\\big|0\\rangle$ to $\\big|1\\rangle$ if $S_T \\geq K_1$ and $S_T \\leq K_2$, and these ancillas are used to control the linear part of the payoff function.\n",
    "\n",
    "The linear part itself is then approximated as follows.\n",
    "We exploit the fact that $\\sin^2(y + \\pi/4) \\approx y + 1/2$ for small $|y|$.\n",
    "Thus, for a given approximation rescaling factor $c_\\text{approx} \\in [0, 1]$ and $x \\in [0, 1]$ we consider\n",
    "\n",
    "$$ \\sin^2( \\pi/2 * c_\\text{approx} * ( x - 1/2 ) + \\pi/4) \\approx \\pi/2 * c_\\text{approx} * ( x - 1/2 ) + 1/2 $$\n",
    "\n",
    "for small $c_\\text{approx}$.\n",
    "\n",
    "We can easily construct an operator that acts as\n",
    "\n",
    "$$\\big|x\\rangle \\big|0\\rangle \\mapsto \\big|x\\rangle \\left( \\cos(a*x+b) \\big|0\\rangle + \\sin(a*x+b) \\big|1\\rangle \\right),$$\n",
    "\n",
    "using controlled Y-rotations.\n",
    "\n",
    "Eventually, we are interested in the probability of measuring $\\big|1\\rangle$ in the last qubit, which corresponds to\n",
    "$\\sin^2(a*x+b)$.\n",
    "Together with the approximation above, this allows to approximate the values of interest.\n",
    "The smaller we choose $c_\\text{approx}$, the better the approximation.\n",
    "However, since we are then estimating a property scaled by $c_\\text{approx}$, the number of evaluation qubits $m$ needs to be adjusted accordingly.\n",
    "\n",
    "For more details on the approximation, we refer to:\n",
    "<a href=\"https://arxiv.org/abs/1806.06893\">Quantum Risk Analysis. Woerner, Egger. 2018.</a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set the strike price (should be within the low and the high value of the uncertainty)\n",
    "strike_price_1 = 1.438\n",
    "strike_price_2 = 2.584\n",
    "\n",
    "# set the approximation scaling for the payoff function\n",
    "rescaling_factor = 0.25\n",
    "\n",
    "# setup piecewise linear objective fcuntion\n",
    "breakpoints = [low, strike_price_1, strike_price_2]\n",
    "slopes = [0, 1, 0]\n",
    "offsets = [0, 0, strike_price_2 - strike_price_1]\n",
    "f_min = 0\n",
    "f_max = strike_price_2 - strike_price_1\n",
    "bull_spread_objective = LinearAmplitudeFunction(\n",
    "    num_uncertainty_qubits,\n",
    "    slopes,\n",
    "    offsets,\n",
    "    domain=(low, high),\n",
    "    image=(f_min, f_max),\n",
    "    breakpoints=breakpoints,\n",
    "    rescaling_factor=rescaling_factor\n",
    ")\n",
    "\n",
    "# construct A operator for QAE for the payoff function by\n",
    "# composing the uncertainty model and the objective\n",
    "bull_spread = bull_spread_objective.compose(uncertainty_model, front=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "tags": [
     "nbsphinx-thumbnail"
    ]
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAE+CAYAAAB1DJw3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd5xU1f3/8ddb0QiiBBtWFkv82VIhMZZIsYMNUFExSjRib19Loqhgwdh7RbElGOy9IB27AVuioqICgiiiIOIiUj6/P84ZGcbZ3Zlldu7c3c/z8ZjH7r1z78x7h2E+c++55xyZGc4551wxVkg6gHPOufTx4uGcc65oXjycc84VzYuHc865onnxcM45VzQvHs4554rmxcNVNEkDJFnW7TNJD0naNKE8q0oaKumrmKdPXH+UpE8kLZI0poZ9O+X8LZnbonL+DTHL5vG1/XnO+j4xU8tyZ3Lp0izpAM4V4Btgj/j7JsCFwEhJW5vZd2XOciywN3AYMB34SNK6wM3ADcADwOw6HqM38HHWchKdrTYH+gN3AXOy1j8FbAdUJ5DJpYgXD5cGi8zslfj7K5KmAs8DXQkf1uW0BfC+mT2UWSFpR2BF4A4ze7uAx3jbzP7XUAGXh5l9CXyZdA5X+fy0lUujCfFnu3ga6QZJ70uqjqeObpS0emZjSffnO5UUT9t8IWmluLyWpLvjKalqSWMkdcjafjJwJPDbrFNOAwiFDOCt7FNZxYrP92DOusyprm3icru4fKCkWyV9I2mapPMlrZCz768kPSFpjqR5kl6TtKukTsATcbNP4uNNjvv85LRVXa9L5rWRdIWkU2Oe2fH03jKnxVzj4cXDpVG7+PNzoAXhW38/YE/gXKALyx6RDAZ2krRxZoUkAYcD/zKzhXH1o8DuwOlAL8L/j9GSNov3dweeBiYSTu1sB9wOHB/v7x3XPVVH/hUlNcu61ef/4WXAPGB/4F/AefH3zN+3BfAisB5wTMz+CLAR8Hr8GwF6xMzda3muul6XjAOBnYG+wN+AvYCL6/G3uTQwM7/5rWJvwABgFuEUazPCufrRwFxgvTzbNwN2ILQjtI3rVgCmAudnbdclbrNNXN4jLnfM2mZVwimcW7PW3QWMz3nOTtmPVcvfktku93ZRvH8M8GBtj00onAbck7Pdm8DQrOV/A9OA5jVk2Ss+Truc9X3i+pZFvi6TgY+AZlnrrgE+T/o95LeGuXmbh0uDNYGFWctTgV5mNgNA0p+B/wN+Qfhgy9gcmGpmSyTdCRwmaYCFT7Y+hCKQaXv4AzDTzMZmdjaz7yQ9CexY4r/nIMIHbcZn9XiM53KW3wXaZi13IRxVza/HY2cr5nUZbWbZV469C6wjaSVbenTnGgkvHi4NvgF2IXwD/hz4LBYAJHUH7iFc7XQ28DXhVM0jwCpZj3En4ZRWZ0n/AXqy9NQNcZ+ZeZ77C2CNUv4xwDu2/A3mc3KWf2DZv3dNYMZyPgcU97rkyyTgZyxb/F0j4MXDpcEiMxtfw30HAK+a2XGZFZI65m5kZpMljSAccWxMOJX176xNZgDr5Hn8NoSCVA7fAyvnrGtdz8f6ivDBv7wq4XVxFcgbzF3aNQcW5KzrXcO2gwlHHMcBj5pZ9jflVwmnWHbKrJDUAugGvFC6uLWaRrgUONtu9XyskcCBklap4f4f4s+a7s+ohNfFVSAvHi7thhOupOonaRdJVxGu+MnnUcK3+98RTmP9yMyGAS8B90k6XNJehCurmgOXN1j6ZT0C/ELS1fFvGcjSzpHFOh9oBYyT1Cs+3hmSjoj3vx9/Hi1pW0m/zPcgFfK6uArkxcOl3a3AlcDJwMNAFXBIvg3NbAHwDPApMCLPJvsRitE1hEt9BXQxs0mlj50331OEdpv9CYWkivB31eex3ic0aM8iXE78SHzcKfH+KYQ2nx6ES3qfyP9IQMKvi6tMiu2OzjV6kpoRPjzvMLNzk87jXJqV/chD0maxZ+zbkhbn6/mbZ5/fS7pT0qTYw/V9Sf1rOZ/r3I8krSzp94TTLGsSjlacc8shiauttiaMSfQKsFKB+/QCNgUuBT4EfkUYHO9XhAZQ52qzPvAa4ZLTo81sWsJ5nEu9sp+2krSCmS2Jvz8IrGVmnerYZy0zm5Wzri/hG2S7eP7WOedcmZT9tFWmcBS5z6w8q9+IP9dfvkTOOeeKleZOgtsBS1h2mIe81lprLWvXrl29nuS7775j1VVXrXvDCpGmvGnKCunKm6askK68acoKy5d3woQJs8xs7Xz3pbJ4KEy+cw7wTzPLN3RC5rRWX4A2bdpwxRVX1Ou55s2bR8uW6ZlULU1505QV0pU3TVkhXXnTlBWWL2/nzp1rbhJIclRG4EFgTJH7rAyMI8zE1rqQfdq3b2/1NXr06Hrvm4Q05U1TVrN05U1TVrN05U1TVrPly0vOCNLZt1QdecQ5GO4hXLG1g5nVNd2nc865BpCq4kHo4bovsKuZTUw6jHPONVWpKR6SzgJOAA40Mx+QzTnnElT24hFH5OwaFzcAVpeUmT7zaTOrljQJGGtmR8Z9DiFMZ3kXMF3SH7Me8iMz+7I86Z1zzkEyAyOuQxhc7QHgj8BWWcuZeQOaEealzsgMS90HeDnn1q3BEzvn3JAh0K4dHbt0gXbtwnIla+C8ZT/yMLPJhFE5a9umXc5yH0LhcM658hsyBPr2herq8OE1ZUpYBuhd0/QxCSpD3tS0eTjnXGL69YPq6mXXVVfD8cfD++/n3ydJ112XP2+/fl48nHOubKZOzb/+m2/goovKm6UQNY1ZWNPfUQ8+GZRzztVlzTXzr6+qgiVLKu9WVZU/b9u2JXtJvHg451xtrr0WZs2CFXI+Llu0gIEDk8lUl4EDQ75sJc7rxcM55/Ixg7PPhlNOgR49YPBgqKrCpPDNftCgymwsh5Br0KAGzettHs45l2vRIjj6aLjjjvDzxhthxRWhTx/GjhlDp06dkk5Yt969oXfvBsvrRx7OOZdt/nzo2TMUjvPOg5tvDoXDLcOPPJxzLmP2bNhnH3jxxXC0cdxxSSeqWF48nHMOYPp02GMP+OADuO8+OOCApBNVNC8ezjn3/vuw227hyOOZZ6BLl6QTVTwvHs65pu2116Br19CuMWYM/O53SSdKBW8wd841Xc89F44yWrUK7RxeOArmxcM51zTdey906wabbRYKx2abJZ0oVbx4OOeanmuvDf0gdtwRxo6FdddNOlHqePFwzjUdZnDWWUt7jT/zTDhl5YrmDebOuaahpl7jrl78yMM51/hVV4cjjTvugP79vdd4CfiRh3OucZs9G/beG156yXuNl5AXD+dc4+W9xhuMFw/nXOPkvcYblBcP51zj473GG1zZG8wlbSbpVklvS1osaUyB+7WSdKek2ZK+kTREUg1zQzrnmqxhw7zXeBkkcbXV1kBX4H3ggyL2ux/oBPwV6AP8Hni0xNmcc2l2772w117ea7wMkjht9YSZPQYg6UFgrbp2kLQdsBvQ0czGxXXTgVcl7WJmIxoysHMuBa69NnT+69QJHn3UO/81sLIfeZjZknrstifwRaZwxMd5Dfgk3ueca6q813gi0tJgvgUwMc/69+J9zrmmyHuNJ0ZmltyTx9NWZtapju2GA9+Z2X456/8FbGJm2+fZpy/QF6BNmzbthw4dWq+M8+bNo2XLlvXaNwlpypumrJCuvGnKCvXLu8L337PVhRey1ksvMfnww5l8+OEgNVDCpZrCa5vRuXPnCWbWIe+dZpbYDXgQGFPAdsOBR/Os/xfwUl37t2/f3upr9OjR9d43CWnKm6asZunKm6asZvXI+/XXZjvsYCaZ3Xhjg2SqSaN/bbMA462Gz9W0nLaaDaydZ33reJ9zrqmYPh123x0+/BDuvx/23z/pRE1SWgZGnEj+to2a2kKcc43RxImw/fYwdWpoGPfCkZi0FI9ngHUl7ZhZIakDsEm8zznX2L32Wpi86fvvQ69xH24kUWU/bSWpBaGTIMAGwOqSMl8fnjazakmTgLFmdiSAmb0s6TngHkmnA0uAS4EXzPt4ONf4DRsGPXtCmzbhd+/8l7gk2jzWAR7IWZdZ3hiYTMiVe71dL+Bq4A7CEdOTwEkNltI5VxmGDIE+fWCbbcKpKp8ytiKUvXiY2WSg1uvpzKxdnnVzgL/Em3OuKbjmGjj1VO81XoHS0ubhnGtKMr3GTz01nK7yXuMVJy2X6jrnmopFi6BvX7jzTjjmGLjhBu81XoH8yMM5Vzkyc43feWeYa/ymm7xwVCg/8nDOVYbsucZvugmOPTbpRK4WXjycc8kZMgT69aPj1KnQrBksWeK9xlPCi4dzLhlDhoS2jerqcPnlwoXws5/BggVJJ3MF8DYP51wy+vULbRzZFiwI613F8+LhnEvG1KnFrXcVxYuHc678hgwJfTnyadu2vFlcvXjxcM6V1zXXwKGHwpZbQvPmy97XogUMHJhMLlcULx7OufLI7jXeowe8/jrcdhtUVWESVFXBoEHQu3fSSV0BvHg45xreokXw17/CJZeEucbvvx9WWSUUismTGTtqFEye7IUjRbx4OOcaVqbX+B13hF7jN9/svcYbAe/n4ZxrONm9xm+8EY47LulErkS8eDjnGsb06bDHHvDBB3DffXDAAUknciXkxcM5V3rvvw+77RaOPJ55xqeMbYS8eDjnSuu116Br19CuMWYM/O53SSdyDcAbzJ1zpfPcc+Eoo1UrePFFLxyNmBcP51xp3HsvdOsGm20WCsdmmyWdyDUgLx7OueV37bWhj8YOO8DYsbDuukkncg3Mi4dzrv7M4Oyz4ZRTQl+OZ5/1ucabCG8wd87Vz6JFYY7xwYPDvBw+ZWyTUvYjD0lbSRopqVrSZ5IukFTnO05SB0nPSfo63kZI2rYcmZ1zOebPh549Q+E47zy45RYvHE1MWYuHpNbACMCAfYELgNOA8+vYb6O4XzPgz/HWDBguqaohMzvncsyeHfpwPPEE3HADnH8+SEmncmVW7tNWxwDNgR5mNpfw4b86MEDSZXFdPt2A1YDuZvYNgKSXgFlAV+Dmho/unOOzz2D33UMnwKFD4cADk07kElLu01Z7AsNyisRQQkHpWMt+KwGLgO+y1s2L6/wrj3Pl8MEHsP32YfTbZ57xwtHElbt4bAFMzF5hZlOB6nhfTR6K21wpaR1J6wBXA7OBBxooq3Mu4z//CZfhVleHXuM775x0IpcwWU1TQTbEk0kLgTPM7Jqc9dOAe8zs7Fr2/Q3wJLBBXDUD2NPM3qph+75AX4A2bdq0Hzp0aL0yz5s3j5YtW9Zr3ySkKW+askK68pYya+vx49nm3HP5oXVr3r7sMuZvuGFJHjdbU31ty2F58nbu3HmCmXXIe6eZle0GLAROybN+GnBxLfutB3wIPAbsEW9PxP3a1vW87du3t/oaPXp0vfdNQprypimrWbrylizrvfearbSS2a9/bfbZZ6V5zDya5GtbJsuTFxhvNXyulvu01WwgXw+i1vG+mpxBaPfY38yeNbNngZ7AYuD0kqd0zsF118Ehh4R2jrFjYb31kk7kKki5i8dEcto24mW4LchpC8mxBfCOmS3MrDCzH4B3gE0bIKdzTZcZ9OsHJ58M3bt7r3GXV7mLxzPA7pJWy1rXC5gPjK1lvynANpJWzqyQ9DNgG2ByA+R0rmlatAiOOgouvjj0Gn/ggTDXuHM5yl08bgEWAA9L2iU2ag8ArrKsy3clTZI0OGu/24H1gUckdZO0F/AooS1kUNnSO9eYzZ8P++8feo2fe673Gne1KmsnQTObLWln4AZCg/ccwiW3A/LkWjFrvwmS9gD6A/+Mq/8L7Go1XG3lnCvCnDmwzz7wwgtw/fVwwglJJ3IVruwDI5rZu0Ctc1KaWbs860YCIxsolnNN12efhbnGJ070XuOuYD6qrnNN2QcfhOFGZs2Cp5+GXXZJOpFLCS8ezjVV48fDnnuGQQ3HjIH27ZNO5FLEJ4NyrikaPhw6d4aWLcOUsV44XJG8eDjX1AwdGuYa32QTeOkl+MUvkk7kUsiLh3NNyfXXh17j223nvcbdcvHi4VxTYAbnnAMnnQT77gvDhsHPf550Kpdi3mDuXGO3aBEceyzcfnvoPX7TTdDM/+u75eNHHs41ZvPnwwEHhMJxzjlw661eOFxJ+LvIucZkyBDo14+OU6fCBhvAqquGvhzXXQcnnph0OteIePFwrrEYMiQMZlhdHeZmnjYtrD/hBC8cruT8tJVzjUW/fmGa2FxPPFH+LK7Rq7N4SBolaYv4+2GS1mz4WM65ok2dWtx655ZDIUcefwIy1/TdiU++5FxlWnvt/Ovbti1vDtckFNLm8SlwgKR5gICN4+95xVFznXPldN99YXBDKfTpyGjRAgYOTC6Xa7QKOfL4B3AS8BZgwL2EuTRyb/+LP51z5XTDDXDwwbDDDuFS3KoqTIKqKhg0CHr3Tjqha4TqPPIws9skPQ78AhgHHA/40YVzSTOD886Diy4Kvcb//W9o3hyOOoqxY8bQqVOnpBO6RqzO4iHpMOApM3tB0vnAY2b2WcNHc87VaNEiOO44uO02OPLIMGWsd/5zZVTIaavsRvLzgA0bLo5zrk7ffx96jd92W7g897bbvHC4sivkHTcbWD/+LkK7h3MuCXPmhFNU48Z5r3GXqEKKxwjgn5Lej8t3Sfqupo3N7A8lSeacW9aMGWGu8ffeC+0bBx2UdCLXhBVSPI4AjgW2AH4HfAJ82ZChnHM5PvwwzDU+cyY89RTsumvSiVwTV8jVVtXAlQCSdgH6mdlb9X1CSVsB1wPbAXOA24HzzWxxAfv2AM4CtgGqgf8APc2sxiMh51JvwoQw17gZjB4Nv/990omcK25gRDPbeHmeTFJrwmmwd4F9CQ3xVxIa7s+pY9+/AjcAlwFnAK2BLvjgjq4xGzECuneHNdeE556DzTdPOpFzQD0+eCVtQvjw3hFYA/gaeB64wsw+rmP3Y4DmQA8zmwsMl7Q6MEDSZXFdvudcC7gaONHMbsu665Fi8zuXGvffD4ceCltsAc8+C+uvX/c+zpVJUaPqSmoPvAn0JJwyuif+7Am8Iel3dTzEnsCwnCIxlFBQOtay34Hx593F5HUutW68MTSIb7ttuLLKC4erMMUOyX4F8AbQzsyOMLOzzOwIYOO4/oo69t8CmJi9wsymEtovtqhlv22B94EjJU2TtFDSq5K2LzK/c5Ut02v8hBNg773DqSqfa9xVIJkV3m0jXqJ7oJk9lee+vYD7zGzVWvZfCJxhZtfkrJ8G3GNmZ9ew3zBge2AucCbwVfzZAfiFmX2RZ5++QF+ANm3atB86dGhhf2SOefPm0bJly3rtm4Q05U1TVihD3sWL2fyaa1j/ySeZseeefHDaadiKK9brofy1bThpygrLl7dz584TzKxD3jvNrOAbMAs4rIb7DgO+qmP/hcApedZPAy6uZb/nCJ0T98hatzqhA+OFdeVu37691dfo0aPrvW8S0pQ3TVnNGjjv/Plm3bubgdlZZ5ktWbJcD+evbcNJU1az5csLjLcaPleLbTB/CrhE0sdm9kJmpaQdCaPv1jVl2WygVZ71reN9te1nwJjMCjObK2kCsFVh0Z2rUN98E3qNjx0L11wDJ5+cdCLn6lRs8fg/4DFgrKSZwExgnXh7GTitjv0nktO2IWkjoAU5bSE53iMMjaKc9QKWFBreuYozY0bow/HOO2EO8kMOSTqRcwUpqsHczL4ysx2BbsBNwIvx555mtqOZfVXHQzwD7C5ptax1vYD5wNha9nsy/uycWSGpFdCeMM+Ic+kzaVKYg2PSpNBr3AuHS5GijjwkrWhmi83sWeDZejzfLYSJpR6WdCmwCTAAuMqyLt+VNAkYa2ZHApjZeEmPAYMl/Z3Q9nImoQ3lxnrkcC5Zr78exqlasgRGjYI/+JBwLl2KvVR3uqTLJG1Znyczs9nAzsCKhPaR8wmd//rnbNosbpPtUOBR4CrgQULh6BIf07n0GDkSOnYMEze9+KIXDpdKxbZ53EK4quo0SeOBwcBQq6FneD4W5jjvUsc27fKsm0cYoPHYYgI7V1EeeCD0Gt9889BrfIMNkk7kXL0U2+YxwMw2AXYldNq7CpghaUgcNNE5V5ObboJevcLAhuPGeeFwqVbsaSsAzGyUmR0GrAucCPw/YJikyZIGSPKxFJzLMIP+/eH442GvvWD4cGjdOulUzi2XehWPLB2AnQiX384mDJD4V2CSpEOX87GdS7/Fi+HYY+GCC+CII+Dhh0Nbh3MpV3TxkFQlqb+kj4CRwHqECaPWN7M/A1XArcDlJU3qXNp8/z0ceCDceiv8/e9w++0+17hrNIq9VHc08CdgOnAncKeZTcnexswWS7oX8G6yrun65hvYbz8YMwauvhpOOSXpRM6VVLFfg2YCXYHhcdyTmrxJGGnXuabn889DH4533oF//Qt69046kXMlV+xMgr0K3G4hMKXODZ1rbCZNCnONf/45PPFEKCLONUL1OgEraUNgc2CV3PvM7OnlDeVcKr3+ehinavHi0Gt8222TTuRcgym2zWM14H5gt8yq+DP7FFb9JiBwLs1GjQptHK1bw7BhYepY5xqxYq+2+gfQltBoLqA70InQ0/wT4I+lDOdcKjz4YDjiaNs2DDfihcM1AcUWj67AQODVuPyZmY0zs76EodrPKGU45yrezTeHy3E7dAi9xjfcMOlEzpVFscWjDfCpmS0GvgPWyLrvaZaeznKucTODAQPguOOgW7fQa3yNNerczbnGotji8SmwVvz9Q2CvrPu2Bb4vRSjnKtrixaFonH8+9OkDjzwCLVokncq5sir2aqvhwC7AI4Sh1O+W1B5YQBim5MrSxnOuwnz/fRgV96GH4G9/g3/8A5Q7waVzjV9BxUNSc0J7x0zgY0ltzOyfkuYB+wPNgRMIw5I417gMGQL9+tFx6lRYeWVYsACuugpOPTXpZM4lps7iIWkTYATQLmv1XEkHmtkjhKMQ5xqnIUOgb1+org7XpS9YEArIOuskncy5RBXS5nEZsIRweW4LYGvgDfwowzUF/fpBdfWy6374Iax3rgkrpHhsB5xjZi+a2fdm9h5wNNBW0noNG8+5hE2dWtx655qIQorHesDHOes+InQSXLfkiZyrFKNH13xf27bly+FcBSr0Ut3aRtB1rvF58MEwqOF66/108qYWLWDgwGRyOVchCi0ewyTNzNyAGXH9yOz18T7n0u2WW5b2Gv/vf+G226CqCpOgqgoGDfJh1l2TV8iluuc3eArnKoFZmC52wIDQa/z++8NRRu/e0Ls3Y8eMoVOnTkmndK4i1Fk8zKykxUPSVsD1hIb4OcDtwPlxyJNC9l8BeA1oD+xtZk+WMp9rohYvhhNPDGNVHX54ONpYaaWkUzlXsco6obKk1oQ+I+8C+wKbEnqlrwCcU+DD/BXw0edc6SxYEHqNP/ggnHkmXHKJ9xp3rg5lLR7AMYTe6D3MbC4wXNLqwABJl8V1NYrFZyDwd8IRi3PLZ+7cMA/H6NFwxRVw2mlJJ3IuFYodGHF57QkMyykSQwkFpWMB+18IvAiMbIBsrqn54gvo1Amefx7uuccLh3NFKPeRxxbAqOwVZjZVUnW874madpT0K+AI4FcNmtA1DR9/DLvtBjNmwOOPh8mcnHMFk1n5unBIWgicYWbX5KyfBtxjZmfXsu9Y4FUzO1NSO8LMhTU2mEvqC/QFaNOmTfuhQ4fWK/O8efNo2bJlvfZNQpryJpW15aRJ/OrMM9Hixfz3H/9g7lZbFbSfv7YNJ01505QVli9v586dJ5hZh7x3mlnZbsBC4JQ866cBF9ey30HA58DqcbkdoePiXoU8b/v27a2+Ro8eXe99k5CmvIlkHT3abLXVzDbayOzdd4vcdXSDRGoIacpqlq68acpqtnx5gfFWw+dquds8ZgOt8qxvHe/7CUkrAZcDlwIrSPo5sHq8e1VJqzVEUNcIPfww7L57mCr2xRdhyy2TTuRcapW7eEwktG38SNJGhNF6J9awz6qES3OvIhSY2cBb8b6hhBF+navdrbfCAQdA+/bwwguw0UZJJ3Iu1crdYP4McIak1czs27iuFzAfGFvDPvOAzjnr1gX+DZxNTgO8c8swgwsvhP79oWtXeOABnzLWuRIod/G4BTgJeFjSpcAmwADgKsu6fFfSJGCsmR1pZouAMdkPEhvMAf5rZq82fGyXSosXw8knw403wmGHwe23e69x50qkrKetzGw2sDOwIuGy3PMJc6H3z9m0WdzGufpZsAAOPjgUjjPOgLvu8sLhXAmV+8gDM3sX6FLHNu3quH8yYT4R535q7lzo3h1GjYLLL4fTT086kXONTtmLh3MN6osvQtvGW2/B3XeH01XOuZLz4uEaj48/DpfiTp8eeo137Zp0IucaLS8ernF4660w89+CBTByJGy3XdKJnGvUyt3Pw7nSGzsWdtoJmjULfTi8cDjX4Lx4uHTL9BrfYAN46SUocJwq59zy8eLh0mvQoNBr/Le/DcOqe69x58rGi4dLn0yv8aOPDkcdI0bAmmsmncq5JsUbzF26LFkCJ50UOv/9+c8weLB3/nMuAX7k4dIju9f46ad7r3HnEuRHHi4dvv029BofOdJ7jTtXAbx4uMo3c2bo8Pfmm95r3LkK4cXDVbZPPglzjU+fDo89Bt26JZ3IOYcXD1fJvNe4cxXLG8xdZRo3znuNO1fBvHi4yvPoo+FU1frre69x5yqUFw9XWW6/HXr2hN/8xucad66CefFwlcEMBg6Eo44KvcZHjvRe485VMG8wd8lbsiTMNX7DDXDooXDHHd75z7kK58XDld+QIdCvHx2nTg2npdZfH155BU47DS67DFbwA2LnKp0XD1deQ4ZA375QXR0moZ86NdwOOgiuuCLpdM65AvlXPFde/fpBdfVP17/8cvmzOOfqrezFQ9JWkkZKqpb0maQLJK1Yxz6/l3SnpElxv/cl9Ze0SrlyuxKZOrW49c65ilTW01aSWgMjgHeBfYFNgSsJReycWnbtFbe9FPgQ+BVwYfzZswEju1Jbd12YMeOn69u2LX8W51y9lbvN4xigOdDDzOYCwyWtDgyQdFlcl88lZjYra3mMpO+BWyVVmUy69VEAABuvSURBVNmUBs7tSmHcOJgzB6RwaW5GixbhMl3nXGqU+7TVnsCwnCIxlFBQOta0U07hyHgj/ly/dPFcg8n0Gq+qgmuugaoqTArLgwZB795JJ3TOFaHcxWMLYGL2CjObClTH+4qxHbAE+Kg00VyDye01ftJJMHkyY0eNgsmTvXA4l0Ky7NMHDf1k0kLgDDO7Jmf9NOAeMzu7wMdZF3gbeNrM+tSwTV+gL0CbNm3aDx06tF6Z582bR8uWLeu1bxIqKq8ZbYcMYZPBg/lq2215p39/ljRv/uPdFZW1AGnKm6askK68acoKy5e3c+fOE8ysQ947zaxsN2AhcEqe9dOAiwt8jJWBccDHQOtC9mnfvr3V1+jRo+u9bxIqJu/ixWYnnmgGZoceavbDDz/ZpGKyFihNedOU1SxdedOU1Wz58gLjrYbP1XI3mM8GWuVZ3zreVytJAu4BtgZ2MLM693EJWLAADj8c7rvPe40710iVu3hMJKdtQ9JGQAty2kJqcA3hEt9dzayQ7V25ffst9OgBI0aEonHGGUkncs41gHIXj2eAMyStZmbfxnW9gPnA2Np2lHQWcAJwoJm90LAxXb3MnBmmiX3jDbjrrnD04ZxrlMp9LuEWYAHwsKRdYqP2AOAqy7p8N/YkH5y1fAhwMeGU1XRJf8y6rV3eP8Hl9cknsOOO8M47Ya5xLxzONWplPfIws9mSdgZuAJ4A5gBXEwpIbq7sIUt2iz/7xFu2vwB3lTapK8rbb4c5OBYsCKertt8+6UTOuQZW9lF1zexdoEsd27TLWe7DT4uGqwTjxsE++8Bqq4UJnHzKWOeaBL8ExtVfptf4euv5XOPONTFePFz9+FzjzjVpXjxccXyuceccPpOgK0b2XON//jMMHuxzjTvXRPmRhyvMggVwyCGhcJx2WujH4YXDuSbLjzxc3bJ7jV9+OZx+etKJnHMJ8+LhajdzJnTtCm++CXffDYcdlnQi51wF8OLhavbJJ+FS3OnTQ6/xbt2STuScqxBePFx+b70Fe+wR2jpGjoTttks6kXOugniDufupsWNhp52gWbPQh8MLh3MuhxcPt6xHHgn9NzbYwHuNO+dq5MXDLXXbbbD//vDb38Lzz3uvcedcjbx4uNBr/KKLoG/fcNQxYoT3GnfO1cobzJu6JUvgpJPgxhu917hzrmB+5NGULVgABx8cCsfpp3uvcedcwfzIo6n69lvo3j1chuu9xp1zRfLi0RTNnAl77hn6cnivcedcPXjxaGqye40//ngYesQ554rkxaMp8V7jzrkS8QbzpsJ7jTvnSsiLR1PgvcadcyXmxaOxGzTIe40750qu7MVD0laSRkqqlvSZpAskrVjAfq0k3SlptqRvJA2R5N2ga2IGF14IRx8d2jm817hzroTK2mAuqTUwAngX2BfYFLiSUMTOqWP3+4HNgb8CS4BLgUeBPzVU3tRavDjMNX7jjeEy3Ntv985/zrmSKveRxzFAc6CHmQ03s1uA84H/k7R6TTtJ2g7YDTjczB4ys0eAQ4EdJe3SIEmHDIF27ejYpQu0axeWK1l23tVXX9pr/M47vXA450qu3MVjT2CYmc3NWjeUUFA61rHfF2Y2LrPCzF4DPon3ldaQIWGQwClTkBlMmRKWK7WA5Oatrg4F4ze/gRW8Wcs5V3rl7uexBTAqe4WZTZVUHe97opb9JuZZ/168r7T69QsfwNmqq6FPH7j44pI/3XL74ANYtGjZdQsXhr+jd+9kMjnnGrVyF4/WwJw862fH++qz3yb5dpDUF+gL0KZNG8aMGVNwyI5Tp6I8623RIr5ce+2CH6dc1n733fx5p05lbBF/d7nNmzevqH+XpKUpb5qyQrrypikrNFzeRtvD3MwGAYMAOnToYJ06dSp857Ztw6mqHKqqYp1KfNO0a5c/b9u2FPV3l9mYMWMqOl+uNOVNU1ZIV940ZYWGy1vuE+KzgVZ51reO95V6v/oZOBBatFh2XYsWYX0lSlte51zqlbt4TCSnjULSRkAL8rdp1LhfVFNbyPLp3Tt0rquqwiSoqgrLldp+kLa8zrnUK3fxeAbYXdJqWet6AfOBsXXst66kHTMrJHUgtHc80xBB6d0bJk9m7KhRMHly5X8Qpy2vcy7Vyl08bgEWAA9L2iU2ag8Arsq+fFfSJEmDM8tm9jLwHHCPpB6S9gOGAC+Y2Yiy/gXOOefKWzzMbDawM7Ai4bLc84Grgf45mzaL22TrRTg6uQO4B5gAdG/IvM455/Ir+9VWZvYu0KWObdrlWTcH+Eu8OeecS5B3P3bOOVc0Lx7OOeeKJjNLOkODk/Ql8NNedIVZC5hVwjgNLU1505QV0pU3TVkhXXnTlBWWL2+VmeUdVqNJFI/lIWm8mXVIOkeh0pQ3TVkhXXnTlBXSlTdNWaHh8vppK+ecc0Xz4uGcc65oXjzqNijpAEVKU940ZYV05U1TVkhX3jRlhQbK620ezjnniuZHHs4554rmxcM551zRvHg455wrmhcP55xzRfPi4ZxzrmiNdg7zxi7OwNgVEPCAmX0laUPgdGBTYDIwyMz+m1xKkPQ34OmkcxRKUnOgmZl9m7VubeAEYCtgCfAmcJOZfZNMSueS55fqRpJEmB+kG7AlsAbhg+Jz4BXgLjP7ILmES0n6AzAcWBVYBHwN7A48DSwG3gG2AdYFdjGz5xOKiqQlgBGmC74XGGpmHyWVpy6SngY+NLOT4/J2hNkqlxDmkBHQHvgB6GJm7ySY9bdAczN7KWvdHsBZLC10bwEDsrepFPH/3N7A7wjvkfGELxoV/aEkaXXCWFFdzOyFpPPAj5m6ACsDT5nZd/FLz/GEGVc/JnyZ/Kxkz1nh/05lEV/kpwkfCp8TPhg2ILyhnyG8+P8PuNDMLkwqZ4ak4YSjxu7Ad4QJtfYjfLjtb2YLJf0MeBRYxcw6J5h1CXAp8EtgV0Lu1wmF5H4zm55UtnwkzQKONLPH4vIrhNd4v8zRiKRWwOPA92a2e4JZXwGeMLOBcfkI4HZgNDCKUOh2Bv4E9Mz8TQllfYnwur4Xl1sTZgdtD8yLm7UkfFHbPfvILwmSjqvl7ubA5cC1wIcAZnZTOXLlI2kzYCSwUVz1CbAb4Qvmz4GPCJ9f84H2ZjatJE9sZk3+Bvyb8Cb4Zda69YFngYfickfCm/yICsj7FbBn1vI6hG+Zu+Vs1w2YlXDWJcAf4u9rAH3jG31RvI2J69ZI+nWNGauBnbKWf8h9XbNe2+8Szjo3OxswCbg+z3a3AG9VyvsgLg8mHDHvkbVuD2A2cHUFvA+WEI7il9Rwy75vccJZ7yccYW4W/4/9M36evQSsFrdZK25za6me1xvMgz2Bv1vWeXkLh3fHAPtJWs/MxgIXAycnlDGX5fk99zCyog4rzexrMxtkZjsDGwKnEQ6zbwFmSHoq0YDB/4DsI7UvCP8hc61JKDRJWpKzXAU8mGe7BwnfPCvJPsAFZvZsZkX8fSDQI7FUSz0GzASOBFY0sxUyN8L7QUCnuC53yuxy2xEYaGaTzOxr4BxCu+cVFo/gzGwWcA3LvreXixePYAXCN4lciwlvklZx+VVg83KFqsV44HRJLSWtAJwNTAeOlbQigKRmwHGED8OKY2afm9m1ZrY9sDFhHvv1E44FcAnwd0lHxNdwIHC5pF0lrSzpZ7Fd4R+Eb3xJeh7onbX8DpBv6O3fE94fleTnhPdxrgmEtrpEmVl34HDgDOA/knbIvjuZVDVqTTjdnpH5t86dw+hjwpe2kvCrrYLhwEWS3jazj+HHc7LXEf5RMg3lLYFKuMKmHyHzbMKpn2pCY9kDwIeSMg3m6xNOBVQ0M5tC+NC+pAKyPCzpRMK3tKuB9wlfHp4lfGgobvo44YMlSWcDL8YvENcTGsrvlrQG4XSgCO+LU4C/JxUyS09JmeI2G8g3ydBahNNxiTOz5yT9ivD6PSXpWcLVjIm2x+Qxk3DUmbEYuJVw1JxtHUqY3RvMgXiJ6zDCUcUUwnnujYEFwMFm9kzc7jLCzFq9ksqaETPvRfgC8JCZzZC0LnAm4RTFFOB2M3s9wZhI6g/cZiW8yqMcJK0J9AL+QPgmLMIH3nvAk2Y2IcF4P5L0G+BmYFuWLW6Z32cTTg9dm0zCIF44kesuMzsiZ7tbga3M7E/lSVaY+H/rMsIptVsJBaWzmY1LNBgg6VHg69zXMs921wNbmNmuJXleLx5BPN1zIPBrYBVC4+O98RyicxVN0paEApJb6F4ys4VJZiuGpKOAj8xsVNJZ8omXbl9N+ILWzSrgEmhJbYAWZvZJHdv9H+HCiZEleV4vHo2PpBXMLN83vYohaRVCo94SYFIlfsDFNo9NyOrzY2ZTk03lXGXwBvMckraW1FPSX+Otp6Stk86VS1IPSY9KelrS3nFdL0mTgYWSpsRvcYmSdGjsf5BZbibpEsJlmm8TGvS/llQJ5+QBkNRe0uOE88PvAS8CLwOfSJou6QJJLRIN2YgoSjpHPpKa5/5bS/pN/Fxon1SuipDk9cmVdAOOILQT5Lu2ezFhuI+/JJ0zZj0w5nqBcElhNXAUoa1mMKFX6b9j7t0TzvoucGzW8pUx77nADoTLDAcQOjCdXQGv7W6Etq7xhEuzBxA6iv4QM59GuKrpTaB1BeTdi9Bv5r34XtgpzzbbknxfhN2IfQ6y1u1H6DC6CFgYX/NuSb+mMVsr4JGYaxFwG7AicHfO58ILwFpJ5y3wb+pZyvdB4n9QJdyAE+Ob5EZCb9y14htlxfj7jsAN8UPl+ArI+x/glqzl3jHblTnb3QmMSDhrNdAxa3kmcHKe7U4HplTAazsBuLuG98hkwtH6KvFD76aEs+4aP8BejO/PCXH5SuIp6bhdJRSPxSzbSbB7/AB+Kf7bn0Y4ultEnk6ZCeS9jjAEyYnAYfELw0PAp7EQrk3oHzYduDnpvAX+TSUtHt7mAUj6mPBhfFkd250JHGNmm5QnWY055gI9zGxEXG5FaCDdxbIaGuPprFvNLLH+E5JmACeY2UNxeQHhaGhMzna7Ao+bWfPyp1wmx3xgHzMbnrO+NaFn/9Zm9p6kw4BLzWy9JHLGTC8QxuH6S9a6IwgffMMJVwp+L2lbQsN5Yp3Z4tVWfzSz1+Ly68B0M9s7Z7ungVXNrGMCMbNzTCZ0vLstLv+WUJz/YmZ3Z213FOGIeeNEgoYMdxS4aRWhY2NJ3gfe5hGsB7xWwHavUQEdmAiXYWa/ATJjA83J2W4eoTNWkh4ndGhcOS6PAA7Os93BhG93SZtJuOIu168Jr3umn88UlnYeTco2wL+yV5jZHYShdP4IjIp9PirRNoRLXnMNIgyUmLS1Wdq/C+IYVoRxorJNIn9/lXI6nHA09Ms6blU1PUB9eCfB4C3gKEnjrIarlGKD3lGERt6kTSGMmjoMwMwWx0sI38vZblNgRpmz5TqL0BP6f5JuB54ALpW0DUs7snUGfksYYTVpg4ALJbUkFLofCD20+wGjbWl/lU2ApK+8+p4wsvIyzGxC7BE9jHBaaECZc9Uk+zTHN+TvsPYdlfGl9hNCER4bl/9EOM22PaGdI2MHkn8ffAi8ZmaH1baRpP2B+0r1pF48gtMIPYjflfQwYfjwzLf4VsAWhHO0G1IZPbYfJmeYATN7Nc92B7HsG73szOxrSX8kfPj+H6GXK8B28fYD4RTLn8zsP8mkXMrMBsZTLH8HzsusJlyAcErWpgsJDepJeptw3v3x3DvM7ONYQJ4G7ipzrpoMk7Qo/t6K8IVhbM42W5D8Fx4I461dK+mXhEJ3IOGL0Lnxi8VbhCOkU4GkR9p+hVDU6pLdiXS5eZtHJGlTQu/sPVg6tHHGp4Qrbi63Cp6LIpektsAcM6uI4R4AJLVj2Y5sH1ll9vFYiXDktgrwcSW9hhmSjiYMUfJbq6Ezq6RVCVcN7WJhUL9ExJEGcn1oZvfmbDcmrq+Ey8xPIpxOXYkwWsMtkg4mtCllBsYcBPwtyfdwvGR4BzO7ro7t1iK02eUW7Po9rxePn4rXdWfaCuaYWdKjpzrnKkQ8hb2WmX2ZdJYkefFoZOIh9etA70o4DaQUTuuqlEzx61ySvHhkiR8a6xNOpczKc/9aQFczu6fs4ZbN0bWWu1clNIr9nTgcu5k9XY5c+ShF07pCuqb4LVQc9+oAM7sg4Rxlnyq1lOIRR/a0uRMIf0fiH6IKoxX3JPx/usvMJkr6NXA+S7/w3GBmw0r2pEl3XKmEG/AzwnDmi+NtIaGndquc7RLvbBVzpGmWs1nAvlnLrxB6RK+Wta4VoeF0WAW8tsMJ07j+nHCu+wZgGqH39kpZ75dnCFdfJf7+LeBvKmnnsHpm2IxwlWDmffkR4UPtY0KB/g9hKPYvgA0r4DV7Cdgya7l1zLgk5pzL0k6OqyWVM2bbnfDl6/P4us4lXME4O+a7Mf6/W0yYTrk0z5v0P1Il3AhX1cwhXIrbATgpvok/BH6RtV2lFI/xhCtS/kK4djv79qv4pj4wsy7hrKmZ1jXmSNMUv20LvB2T9PuWhKZKXY68qZk2lzDCwAOEGQ8hXEQxGxics90/gVdK9rxJ/yNVwo1wae4JOevWBcYBXwLbxXWVUjxEmPd7JmHIhI2z7msV3/g/GeMooayvAf2zlj8FDsqz3WHAlxWQ96ucD4i14+u5a852XSugeGSOMuu6VcIR6GfAgVnLVTFXj5zt/gJ8UAHvg9zi8SVwSp7tEh9Wh3Ap8S5Zy61j/i452+1GuACoJM/r/TyCjcjp/Gdmn0vamVCtR0jqTWVcf46Fd8IgSfcDFwH/jRO9XJRssrwuAYZI+hS4h6XTun5FOFWV6SRYCdO6wtIpfl8gHDVlT/E7ykKHzEqZ4vdbYBRwex3b7Ui4DD1JiUyVWkKVPG3ufJbtLJr5PXeonxaEjqUl4cUj+Az4BeFI40cWrt0+SNK1hMPCRBvKc5nZHOAESYMI155/CFxKBc2xbOma1hXSNcXva4R2uadq2yjOnZK0RKZKXU5pmTb3ReA8SR/GLFcQRrP+Wxw149s4/t2ZhGJXEn61FT8OLLaJmXWqZZuzCN+azRIcYK42kg4iTJW5IWEAtMSnyMxQSqZ1hVRN8Xsu0NfMcju15m63E3C+mXUuT7K8GRKZKrW+lKJpcyVtRhhKJ/M+mEw4mn+Q0GN/CtCO8GWos5m9WZLn9eLx42VuvYB/WC3Tzko6hHDu+y81bZO0eEplVWCemS1OOo9zkNxUqQ1NFTJtbuzftQPhCsGRZjY/dnb+K0u/8NxrZtNK9pxePJxzzhWrEkavdA1E0m2SBiedoxBpygrpy+tcqXmDeREk3QasYGZHJp2lQJ1JzxeENGWFFOWVNIJwlmHnpLPUJU1ZIV15S53Vi0dxUvOBAWBmmyWdoVBpygqpyyvS875NU1ZIV96SZvU2j0YsXqK5jpklPVlNndKUFdKX17lSS0vFrAiSVolzZKRFN8KMaGmQpqyQorySVkrL+zZNWSFdeUud1YtHcVLzgeGaBknHS/pI0nxJb0n6c57NfkcFvG/TlBXSlTeJrN7mkUKSCr2mPF+P2LJKU1ZIV97YKfR6whS5bxCmIr1L0r7AoWZWsqEolleaskK68iaV1YsH6frAiHYiDPPxbh3bVcKwFGnKCunKezpwhZn9OG5VHI9tCDBa0l5m9lVi6ZaVpqyQrryJZPUGc0DSIgr7wNgA2Dbp4UkkvQVMNLNedWy3P3BfknnTlDXmSE1eSd8Ce5vZmJz17QjzjaxIGH9rbeAlz1q4NOVNKqu3eQTvAP8zswNquwFXJR00egX4YwHbZQ88mJQ0ZYV05f2GMDDfMsxsMuHUxSzgZeD35Y2VV5qyQrryJpLVjzz4cXCzPcysqo7tehLmtE606EraFNjazB6vY7vmhMtJc4e9Lps0ZY05UpNX0mPAt2Z2aA33NycMjrcnCQ/omaasMU9q8iaV1YsH6frAcC5D0gHAqcBeNQ3oKWlF4GbCgJ4blzNfTo7UZI1ZUpM3qaxePJxzzhXN2zycc84VzYuHc865onnxcE2OpD6SJkj6VtJsSW9IapAr6SRtLmmApJ8XsO0ASZZ1+0zSQ7FNrq5975KUb45t5xqEFw/XpChMJ3w7MAzoARwGPAbs00BPuTnQH6izeETfANvF2+nAb4CRklatY78LgT71zOhc0byHuWtqTgBuNbOzs9Y9Ien8pALlWGRmr8TfX5E0FXge6Ao8kLuxpOZmNt/MPipnSOf8yMM1NT8HPs9daVmXHUpqF08bHSLpn/H01kxJ/XP3k9RF0quSvpf0haSbFOaTRlIn4Im46SfxMScXmXdC/NkuPuZkSVdKOlfSNGBuXP+T01aSqiT9W9IsSdWS3pZ0SNb9q0i6TNKnkhbEAfW6FpnPNVF+5OGamteBE+M3+ifrGPPncuBJYH/CmFf9Jc0ysxsBJG0NPAsMB3oCGwGXAJsQhoN4nTjuEOEU2QxgQZF528Wf2QXvEMKoCMdRw/9hSesQehVXxwyfAtvEjBkPAn8gnFb7CDgQeFxSBzN7s8icronx4uGamuOBR4G7AJP0HvAQYWC5uTnbvmNmR8ffh8UP5LMl3WxmS4BzgSnAPma2GEDS18B9krYzs5clvR/3fyMOF1EnSZn/l5sANwHfAiNyNturjtFSTwVaAe3NbEZcNzLrOXYmTDHQyczGxtXPSdoc6AccUEhW13T5aSvXpJjZ28CWhAbymwjjU50LjM+cbsrySM7yw8D6wIZx+Q/AI5nCET0ELAJ2rGfENYGF8fY+oYD0yioAACMLGGa7C/Bszn7ZdiEczbwoqVnmRigwHeqZ3TUhfuThmhwzW0Boi3gCQNKRhCuwjgSuzdp0Zs6umeX1gKnx5xc5j71Y0lfAGvWM9w3hg90IH+6f2U+HgfjiJ3v91JrAf2q5fy1gXUKRyrU4zzrnluHFwzV5ZjZY0mXAFjl3rVPD8oysn8tsE8cQWhPIO8ZQARaZWV39NQoZU+grQnGrydfAdGC/QoM5l81PW7kmJbZb5K5bm9A+kPuNvnvOcqbRe1pcfhXoHgtG9jbNgBfi8g/xZ7knjxoJ7C6pTS33rwvMM7PxubfyxXRp5Ucerqn5bxzC+jnCaagqwtVI1cDdOdtuHYfrf4hwtdWRwMmxsRzgIsK0n49KupnQFnIpMMzMXo7bZBrMj5Y0FKg2s/82zJ+2jKsJHSCflzSQcLXVlsCqZnYZ4QqxYcBwSZcSrt5andApcRUzO6sMGV2KefFwTc0FwL7AdYR2ic+BlwiN0p/kbHsmsBeheHxP6MV9Q+ZOM3tH0p7AxYTG9LmEeaTPzNpmiqTTgZOAEwlHLe0a4g/LZmZfStoBuAy4BvgZ8CHwj3i/SeoBnA2cArQlnMp6kzAftnO18iHZncsRp+/8hDC155PJpnGuMnmbh3POuaJ58XDOOVc0P23lnHOuaH7k4ZxzrmhePJxzzhXNi4dzzrmiefFwzjlXNC8ezjnnivb/AYfamBgqivfUAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot exact payoff function (evaluated on the grid of the uncertainty model)\n",
    "x = uncertainty_model.values\n",
    "y = np.minimum(np.maximum(0, x - strike_price_1), strike_price_2 - strike_price_1)\n",
    "plt.plot(x, y, 'ro-')\n",
    "plt.grid()\n",
    "plt.title('Payoff Function', size=15)\n",
    "plt.xlabel('Spot Price', size=15)\n",
    "plt.ylabel('Payoff', size=15)\n",
    "plt.xticks(x, size=15, rotation=90)\n",
    "plt.yticks(size=15)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "exact expected value:\t0.5695\n",
      "exact delta value:   \t0.9291\n"
     ]
    }
   ],
   "source": [
    "# evaluate exact expected value (normalized to the [0, 1] interval)\n",
    "exact_value = np.dot(uncertainty_model.probabilities, y)\n",
    "exact_delta = sum(uncertainty_model.probabilities[np.logical_and(x >= strike_price_1, x <= strike_price_2)])\n",
    "print('exact expected value:\\t%.4f' % exact_value)\n",
    "print('exact delta value:   \\t%.4f' % exact_delta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Evaluate Expected Payoff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set target precision and confidence level\n",
    "epsilon = 0.01\n",
    "alpha = 0.05\n",
    "\n",
    "# construct amplitude estimation \n",
    "ae = IterativeAmplitudeEstimation(epsilon=epsilon, alpha=alpha, \n",
    "                                  state_preparation=bull_spread,\n",
    "                                  objective_qubits=[num_uncertainty_qubits],\n",
    "                                  post_processing=bull_spread_objective.post_processing)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "result = ae.run(quantum_instance=Aer.get_backend('qasm_simulator'), shots=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Exact value:    \t0.5695\n",
      "Estimated value:\t0.5694\n",
      "Confidence interval: \t[0.5639, 0.5750]\n"
     ]
    }
   ],
   "source": [
    "conf_int = np.array(result['confidence_interval'])\n",
    "print('Exact value:    \\t%.4f' % exact_value)\n",
    "print('Estimated value:\\t%.4f' % result['estimation'])\n",
    "print('Confidence interval: \\t[%.4f, %.4f]' % tuple(conf_int))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Evaluate Delta\n",
    "\n",
    "The Delta is a bit simpler to evaluate than the expected payoff.\n",
    "Similarly to the expected payoff, we use comparator circuits and ancilla qubits to identify the cases where $K_1 \\leq S_T \\leq K_2$.\n",
    "However, since we are only interested in the probability of this condition being true, we can directly use an ancilla qubit as the objective qubit in amplitude estimation without any further approximation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# setup piecewise linear objective fcuntion\n",
    "breakpoints = [low, strike_price_1, strike_price_2]\n",
    "slopes = [0, 0, 0]\n",
    "offsets = [0, 1, 0]\n",
    "f_min = 0\n",
    "f_max = 1\n",
    "\n",
    "bull_spread_delta_objective = LinearAmplitudeFunction(\n",
    "    num_uncertainty_qubits, \n",
    "    slopes,\n",
    "    offsets,\n",
    "    domain=(low, high),\n",
    "    image=(f_min, f_max),\n",
    "    breakpoints=breakpoints, \n",
    ")  # no approximation necessary, hence no rescaling factor\n",
    "\n",
    "# construct the A operator by stacking the uncertainty model and payoff function together\n",
    "bull_spread_delta = bull_spread_delta_objective.compose(uncertainty_model, front=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set target precision and confidence level\n",
    "epsilon = 0.01\n",
    "alpha = 0.05\n",
    "\n",
    "# construct amplitude estimation \n",
    "ae_delta = IterativeAmplitudeEstimation(epsilon=epsilon, alpha=alpha, \n",
    "                                        state_preparation=bull_spread_delta,\n",
    "                                        objective_qubits=[num_uncertainty_qubits])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "result_delta = ae_delta.run(quantum_instance=Aer.get_backend('qasm_simulator'), shots=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Exact delta:    \t0.9291\n",
      "Estimated value:\t0.9290\n",
      "Confidence interval: \t[0.9264, 0.9316]\n"
     ]
    }
   ],
   "source": [
    "conf_int = np.array(result_delta['confidence_interval'])\n",
    "print('Exact delta:    \\t%.4f' % exact_delta)\n",
    "print('Estimated value:\\t%.4f' % result_delta['estimation'])\n",
    "print('Confidence interval: \\t[%.4f, %.4f]' % tuple(conf_int))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-08-22T01:55:52.763931Z",
     "start_time": "2019-08-22T01:55:52.753702Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<h3>Version Information</h3><table><tr><th>Qiskit Software</th><th>Version</th></tr><tr><td>Qiskit</td><td>None</td></tr><tr><td>Terra</td><td>0.16.0.dev0+28d8c6a</td></tr><tr><td>Aer</td><td>0.6.1</td></tr><tr><td>Ignis</td><td>0.5.0.dev0+470d8cc</td></tr><tr><td>Aqua</td><td>0.8.0.dev0+ce81016</td></tr><tr><td>IBM Q Provider</td><td>0.8.0</td></tr><tr><th>System information</th></tr><tr><td>Python</td><td>3.7.7 (default, May  6 2020, 04:59:01) \n",
       "[Clang 4.0.1 (tags/RELEASE_401/final)]</td></tr><tr><td>OS</td><td>Darwin</td></tr><tr><td>CPUs</td><td>2</td></tr><tr><td>Memory (Gb)</td><td>16.0</td></tr><tr><td colspan='2'>Fri Oct 16 11:03:18 2020 CEST</td></tr></table>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div style='width: 100%; background-color:#d5d9e0;padding-left: 10px; padding-bottom: 10px; padding-right: 10px; padding-top: 5px'><h3>This code is a part of Qiskit</h3><p>&copy; Copyright IBM 2017, 2020.</p><p>This code is licensed under the Apache License, Version 2.0. You may<br>obtain a copy of this license in the LICENSE.txt file in the root directory<br> of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.<p>Any modifications or derivative works of this code must retain this<br>copyright notice, and modified files need to carry a notice indicating<br>that they have been altered from the originals.</p></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import qiskit.tools.jupyter\n",
    "%qiskit_version_table\n",
    "%qiskit_copyright"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
